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Abstract 



In the first part, the stabihty of two-dimensional parallel flow is discussed. A more restrictively 
general stability criterion for inviscid parallel flow is obtained analytically. First, a sufficient criterion 
for stability is found as either — /xi < jjzrfj- < or < in the flow, where Ug is the velocity 

at the inflection point, and fii is the eigenvalue of Poincare's problem. Second, this criterion is gen- 
eralized to barotropic geophysical flows in the (3 plane. Based on the stability criteria, the necessary 
condition for wave-mean flow interaction is also obtained. 

Then, the general stability criteria of two-dimensional inviscid rotating flow with angular velocity 
J7(r) are obtained analytically. First, a necessary instability criterion for centrifugal flows is derived 
as ^'(O — < (or ^' / [VL — ^s) < ) somewhere in the flow field, where ^' is the vortictiy of profile 
and is the angular velocity at the inflection point ^' = . Second, a criterion for stability is found 
as —{ill + l/r2) < /(r) = qtq7 < , where ni is an eigenvalue. The new criteria are the analogues 
of the criteria for parallel flows, which are special cases of Arnol'd's nonlinear criteria. Specifically, 
Pedley's cirterion is proved to be an special case of Rayleigh's criterion. Moreover, the criteria for 
parallel flows can also be derived from those for the rotating flows. The analogy between rotation and 
stratification in inviscid flow is also addressed. These results extend the previous theorems and would 
intrigue future research on the mechanism of hydrodynamic instability. 

Besides, the essence of shear instability is fully revealed within the linear context. The mechanism 
of shear instability is explored by combining the mechanisms of both the Kelvin-Helmholtz instability 
(K-H instability) and resonance of waves. The shear instability requires both a concentrated vortex 
(with speed of Ug ) in the flow and resonant waves to interact with the concentrated vortex. Physically, 
the standing waves (with phase speed = Us ) can interact with the concentrated vortex, so they 
can trigger instability via K-H instability in the flows. While the travelling waves (with Cr ^ Ug ) 
have no interaction with the concentrated vortex, so that they can not trigger instability in the flows. 
The resonance of waves are totally within the linear context. 

In consequence, the above criteria would be helpful for understanding the wave-mean flow in- 
teraction, especially the Rossby wave-mean flow interaction in barotropic flows. According to the 
stable criteria, the necessary condition for wave-mean flow interaction can be obtained. And why the 
disturbed waves can't take energy from the mean flow in the stable flow is revealed. If the flow is 
stable, there is no wave-mean flow interaction at all. This explains why the disturbed waves can't 
take energy from the mean flow in the stable flow. 

In the second part, we report the numerical simulations of the partial-penetrating flow in hor- 
izontal convection within a squire cavity tank at high Rayleigh numbers 10^ < Ra < lO^'^ . The 
partial-penetrating flow was first reported in the experiment by Wang and Huang (2005), which is 
though of an important material to understand ocean circulation energy budget. The fast established 
but slowly steadied flow is simulated, where a shallow and closed circulation cell is obtained numer- 
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ically as partial-penetrating flow for the first time, which is consistent with the experiment. As the 
partial-penetrating flow is shallow, it is seldom affected by the bottom boundary. The depth of partial- 
penetrating circulation satisfies minus 1/5-power law of Rayleigh number. The larger the Rayleigh 
number is, the shallow the partial-penetrating flow is. An objective definition of partial-penetrating 
is given based on this power law. Then, further investigation points out that the Prandtl number 
governs the partial-penetrating flows. As Pr > 6 , there are always partial-penetrating flows. While 
Pr < 4 , the flows tend to be full-penetrating. 

Then, the horizontal convection at high Rayleigh number (Ra) in a rectangle cavity with aspect 
ratio of 1 : 10 is numerically simulated. According to the results within the regime of 10"^ < Ra < 
10^^ , three continues regimes are obtained: linear regime (10^ < Ra < 10^), transition regime 
( 10^ < Ra < 10^ ) and 1/5-power law regime ( 10^ < Ra < 10^^ ). For the flow strength, a 1/3-power 
law of Ra is fitted when Ra is not high enough (10^ < Ra < 10*). However, a 1/5-power law is 
obtained as Ra is high enough ( 10* < Ra < 10^^ ). The 1/5-power law confirms Rossby's analysis 
and implies that 1/3-power law of Ra for Nusselt number by Siggers et al. is over estimation. 

Finally, the critical Rayleigh number for unset of the horizontal convection is also addressed. The 
flow is found to be unsteady at high Rayleigh numbers. There is a Hopf bifurcation of Ra from steady 
solutions to periodic solutions, and the critical Rayleigh number Rac is obtained as Rac = 5.5377x 10* 
for the middle plume forcing at Pr = 1 , which is much larger than the formerly obtained value. 
Besides, the unstable perturbations are always generated from the central jet, which implies that the 
onset of instability is due to velocity shear (shear instability) other than thermally dynamics (thermal 
instability). Finally, Paparella and Young's (2000) second hypotheses about the destabilization of the 
flow is numerically proved, i.e. the middle plume forcing can lead to a destabilization of the flow. 

The report was supported by the National Foundation of Natural Science (No. 40705027) and 
the National Science Foundation for Post-doctoral Scientists of China (No. 20070410213). 
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Chapter 1 1 

General Stability Criteria 



§1.1 Two-Dimensional Problem 

§ 1.1.1 Introduction 

The instability due to shear in the flow is one of the fundamental and the most attracting problems 
in many fields, such as fluid dynamics, astrophysical fluid dynamics, oceanography, meteorology, etc. 
More generally, shear instability is also referred as barotropic instability in geophysical flows, where 
the gravitational and buoyancy effects are ignored. Shear instability has been intensively investigated, 
which is to the greatly helpful understanding of other instability mechanisms in complex shear flows. 
Rayleigh investigated the growth of linear disturbances by means of normal mode expansion, which 
leads to Rayleigh's equation [1]. Using this equation, Rayleigh first proved a necessary criterion for 
instability, i.e., Infiection Point Theorem, which is also called Rayleigh-Kuo theorem [e.g. 2] for Kuo's 
generalization to barotropic geophysical flows in the P plane [3]. Then, Fjortoft found a stronger 
necessary criterion for instability [4]. Besides, Tollmien gave a heuristic result that the criteria are 
also sufficient for instability if the velocity profiles are the symmetrical or monotonic [5] . These criteria 
are well known and have been widely used in various applications [e.g. 2, 6, 7]. 

On the other hand, Arnol'd considered the shear instability in a totally different way [8, 9, 10]. He 
studied the conservation law of the inviscid flow via Euler's equations and found two nonlinear stability 
theorems by means of variational principle. Arnol'd's first stability theorem corresponds to Fj0rtoft's 
criterion [6, 11]. However, Arnol'd's second nonlinear stability theorem, has no such corresponding 
linear criterion. Though Arnol'd's second nonlinear theorem is more useful in the geophysical fiows 
[11], is seldom known by the scientists in other fields. Dowling suggested that Arnol'd's idea should 
need to be added to the general fluid-dynamics curriculum [11]. Yet his suggestion has not been 
followed until now [e.g. 2, 6, 7, 12], since the proofs of Arnol'd's theorems are very advanced and 
complex in mathematics. 

The aim of this section is to flnd the elementary proofs for Arnol'd's theorems, which could be 
used to teach undergraduate students. As the variational method is too advanced and complex for 
undergraduate students. The new proofs are obtained in a totally different way, where the linear 
stability problem is considered by using Rayleigh's equation. 
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§ 1.1.2 Stable criterion 

To find the criteria, Rayleigh's equation for an inviscid parallel flow is employed [1, 2, 6, 7, 13], 
which is the vorticity equation of the disturbance [6, 7]. For a parallel flow with mean velocity U{y) 
in Fig. 1-1, the vorticity is conserved along pathlines. The amplitude of disturbed flow streamfunction, 
namely (p , satisfies 

- fcV) - j^<P = 0, (1-1) 
U — c 

where k is the nonnegative real wavenumber and c = Cr + ici is the complex phase speed and double 
prime " denotes (^^ jdy^ . The real part of complex phase speed is the wave phase speed, and 
uji = kci is the growth rate of the wave. This equation is to be solved subject to homogeneous 
boundary conditions 

(P = Oaty = a,b. (1-2) 

There are three main categories of boundaries: (i) open channels with both a and b being finite, (ii) 
boundary layers with either a or b being infinite, and (iii) free shear flows with both a and b being 
infinite. 

It is obvious that the criterion for stability is = (cj = 0), for that the complex conjugate 
quantities cj)* and c* are also a physical solution of Eq.( 1-1) and Eq.(l-2). 

Multiplying Eq.(l-l) by the complex conjugate (f>* and integrating over the domain a < y < b , 
we get the following equations 



"' + k^n') + K§^Uf]dy = 0, (1-3) 



and 



C^ I „J l0f ^y = 0- (1-4) 



U-c 



Rayleigh used only Eq.(l-4) to prove his theorem, i.e., a necessary condition for instability is U"{ys) = 
, where ys is the inflection point and Ug = U{ys) is the velocity at yg . Fj0rtoft noted that Eq.(l-3) 
should also be satisfied, then he obtained his necessary criterion. To find a more restrictive criterion, 
we shall investigate the conditions for q = . Unlike the former investigations, we consider this 




a 



Figure 1-1: Sketch of parallel flow. 
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1.1 Two-Dimcnsional Problem 



problem in a totally different way: if the velocity profile is stable ( Cj = ), then the hypothesis Cj 7^ 
should result in contradictions in some cases. Following this, some more restrictive criteria can be 
obtained. 

To find a stronger criterion, we need to estimate the ratio of ||</>'|p(iy to ||(/>|p(iy . This is 
known as Poincare's problem: 



dy = fi Ufdy, (1-5) 



where the eigenvalue /i is positive definite for any (/> 7^ . The smallest eigenvalue value, namely /ii , 
can be estimated as ni > (^r^)^ , like what Tollmien have done [5] . 

Then using Poincare's relation Eq.(l-5), a new stability criterion may be found: the parallel flow 
is stable if — /ii < jf^rjj- < everywhere. 

To get this criterion, we introduce an auxiliary function /(y) = jj^jj , where f{y) is finite at 
the inflection point. We will prove the criterion by two steps. First, we will prove proposition 1: if 
the velocity profile is subject to —fii < f{y) < , then Cr ^ Us ■ 

Proof: Since —fii < f{y) < , then 

^ ^" - u"iu-Us) , U"iU-Us) 

Substitution of = Us and Eq.(l-6) into Eq.(l-3) results in 



.,,,2 , .2„^„2 , U"{U-Us) 

\U -c 



k^M?+ ;:,^' Uf]dy>Q. (1-7) 



This contradicts Eq.(l-3). So proposition 1 is proved. 

Then, we will prove proposition 2: if — /^i < f{y) < and Cr ^Us , there must be = . 

Proof: If 7^ , then multiplying Eq.(l-4) by (c^ — Ut)/ci , where the arbitrary real constant 
Ut does not depend on y , and adding the result to Eq.(l-3), yields 



"'' + k'Uf) + K^^Uf]dy = ^. (1-8) 



But the above Eq.(l-8) does not hold for some special Ut ■ For example, if Ut = 2cr — Ug , then there 
is {U - Us){U - Ut) < \\U - cf , and 

U"{U-Ut) „, ^ {U-Us){U-Ut) ^ 
\\U — c||^ \\U — c|p 

This yields 



{\Wf + r + > o> (1-10) 



which also contradicts Eq.(l-8). So proposition 2 is also proved. 

Using "proposition 1: if —^i < f{y) < then Cr 7^ Us" and "proposition 2: if —fii < f{y) < 
and Cr 7^ Us then Cj = " , we find a stability criterion. If the velocity profile satisfies —fii < < 

everywhere, the parallel fiow is stable. Moreover, the above proof is still valid for < /(y) , which is 
equivalent to Fj0rtoft's criterion. Thus we have the following theorem. 
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Theorem 1: If the velocity profile satisfies either —fii < < or < , the flow is 

stable. 

This criterion covers Rayleigh's and Fj0rtoft's criteria. And the proofs here are very simple and 
easy to understand comparing to Arnol'd's proofs. As mentioned above, we have investigated the 
stable criterion via Rayleigh's equation, while Arnol'd [14] considered the hydrodynamic stability in 
a totally different way. Is there any relationship between these proofs? Two points are outlined here. 
First, this criterion is essentially the same as Arnol'd's second stability theorem and is more restrictive 
than Fj0rtoft's criterion. Second, the proofs here are similar to Arnol'd's variational principle method. 
For the arbitrary real number Ut , which is like a Lagrange multiplier in variational principle method, 
plays a key role in the proofs. 

§ 1.1.3 Discussion 

One may note that the criterion above is something different from Fj0rtoft's criterion. Why are 
the functions of U" / {U — Us) used in Arnol'd's theorems and present theorems, unlike U"(U — Ug) in 
Fj0rtoft's theorem? This is due to the property of Rayleigh's equation. It can be seen from Eq.(l-l) 
that the stability of profile U{y) is not only Galilean invariant, but also independent from the the 
magnitude of U{y) due to the linearity. So the stability of U{y) is the same as that of AU{y) + B , 
where A and B are any arbitrary nonzero real numbers. As the value of U" (U — Us) is only Galilean 
invariant but not magnitude free, it satisfies only part of the Rayleigh's equation's properties. On the 
other hand the value of U"/(U — Us) satisfies both conditions, this is the reason why the criteria in 
both Arnol'd's theorems and present theorems are the functions of U" /(U — Ug) ■ Since the stability 
of inviscid parallel flow depends only on the velocity profile's geometry shape, namely f{y) , and the 
magnitude of the velocity profile can be free, then the instability of inviscid parallel flow could be called 
"geometry shape instability" of the velocity profile. This distinguishes from the viscous instability 
associated with the magnitude of the velocity profile. 

Moreover, the above theorem is essentially associated with vorticity distribution in the flow field. 
As known from Fj0rtoft's criterion, the necessary condition for instability is that the base vorticity 
profile ^ = —U' has a local maximum. Note that U"/{U — Ug) ~ Cs/^s near the inflection point, 
where is the vorticity at the inflection point, that means that the base vorticity ^ must be convex 
enough near the local maximum for instability, i.e., the vorticity should be concentrated somewhere 
in the flow for instability. Otherwise, the flow is stable if the vorticity distribution is smooth enough 
near the inflection point at ys ■ A simple example can be obtained by following Tollmien's way [5]. 
Fig. 1-2 depicts three vorticity profiles within the interval —1 < y < 1 , which have local maximal at 
y = . Profile 2 ( f7 = — 2sin(7ry/2)/7r ) is neutrally stable, while profile 1 {U = — sm{y) ) and profile 
3 ( [/ = — sin(2y)/2 ) are stable and unstable, respectively. 

To show the advantage of the criteria obtained above, we consider the stability of velocity profile 
U = tanh(ay) within the interval — 1 < y < 1 , where a is a constant. This velocity profile is an 
classical model of mixing layer, and has been investigated by many researchers (see [2, 7, 13] and 
references therein). Since U"{U -Us) = -2a2 tanh^ (ay) /cosh^ (ay) < for -1 < y < 1 , it might be 
unstable for any a according to both Rayleigh's and Fj0rtoft's criteria. But it can be derived from 
Theorem 1 that the flow is stable for < vr^/S. For example, we choose qi = 1.1 and 02 = 1.3 
for velocity profiles Ui{y) and U2{y) ■ The growth rate of the profiles can be obtained by Chebyshev 
spectral collocation method [13] with 100 collocation points, as shown in Fig. 1-3. It is obvious that 
Cj = for Ui and Cj > for U2 , which agrees well with the criteria obtained above. This is also 
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Figure 1-2: Vorticity profiles within the interval —1 < y < 1 . Profile 2 ( ,^ = cos(7ry/2) , dashed) is 
neutrally stable, while profile 1 ( ^ = cos(y) , solid) and profile 3 ( ^ = cos(2y) , dash doted) are stable 
and unstable, respectively. 



a counterexample that Fj0rtoft's criterion is not sufficient for instability. So this new criterion for 
stability is more useful in real applications. 

Then, recalling the proof of theorem 1, we will find that the following Rayleigh's quotient /(/) 
plays a key role in determining the stability of parallel flows. 



mm 



I' 

J a 



(1-11) 



Note that the proof of theorem 1 is still valid in the case of /(/) > . We have such a result: parallel 
flows are stable if /(/) > . Although this criterion is more restrictive than that in theorem 1, it is 
inconvenient for the real applications due to the unknown Rayleigh's quotient /(/) . Theorem 1 is 
more convenient for the real applications in different research fields. 

The present stability criteria give a affirmative answer to the questions at the beginning, i.e., 
there are some stable flows if U"{U — Us) < ■ Based on the former criteria, the velocity profiles can 
be categorized as follows: (i) without inflection point (Reyleigh's criterion), (ii) /(y) > (Fj0rtoft's 



g- 10- 




0.25 0.5 

k 



9.0X1 



3 6.0X10" 



3.0X10" 




Figure 1-3: Growth rate uJi as an function of wavenumber k , (a) for Ui = tanh(l.ly) , (b) for 
U2 = tanh(1.3y) , both within the interval —1 < y < 1 . 
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U">Oor U"<0 


u;'=o 






A stable 


B f{y)>0 stable 


C f(y)>-n, stable 


D f(y)<-n, unstable 



Figure 1-4: Diagram of stable velocity profiles. 



criterion), and (iii) —^ii < f{y) < (present criterion), which are shown in Fig. 1-4. Then the parallel 
flow might be unstable only for f{y) < —fj,i and /(y) changing sign within the interval. However, if 
f{y) changes sign somewhere within the interval [a, b] , then the parallel flow is stable. For that /(y) 
changing sign implies U'J' = but U'J" 7^ , so U" does not change sign near the inflection point. 
Thus Cj must vanish in Eq.(l-4), i.e., the parallel flow is stable for /(y) changing sign within the 
interval. In this way, the parallel flow might be unstable only for /(y) < fJ-i somewhere, which will 
intrigue further studies on this problem. In fact, there are still stable flows if /ii < /(y) is violated. 

Finally, the stable criterion for the parallel inviscid flows can be applied to the barotropic geo- 
physical flows in a differentially rotating system. Considering the barotropic plane flows in a rotating 
frame, which are the approximations of barotropic geophysical flows [2, 3, 11], Eq.(l-l) changes to 

U" -p 



u 



0, 



(1-12) 



where /? is the gradient of the Coriolis parameter with respect to y. Eq.(l-12) is a generalized 
Rayleigh's equation, and there is a generalized stability criterion for these flows. 

Theorem 2: The flow is stable, if the velocity profile satisfies either —^1 < < or < 

^Jif , where Ug is the velocity at the point U"{ys) = P ■ 

The above criteria would be helpful for understanding the wave-mean flow interaction, especially 
the Rossby wave-mean flow interaction in barotropic flows. According to the stable criteria, the 
necessary condition for wave-mean flow interaction can be obtained. And why the disturbed waves 
can't take energy from the mean flow in the stable flow is revealed. If the flow is stable, there is no 
wave-mean flow interaction at all. 

First, when the velocity profile has no inflection point, then the speeds of barotropic waves 
Cr — U 7^ . According to Eq.(l-3), U — Cr > holds for U" < , i.e., the barotropic waves are 
always west-propagation relative to the mean flow U if the vortex gradient is positive. This extends 
the west-propagation theory of Rossby waves. And U — Cr < holds for U" > , i.e., the barotropic 
waves are always east-propagation relative to the mean flow U if the vortex gradient is negative. 
Second, when the velocity profile has an inflection point U'J = , the speed of most favorite wave 
which might have interaction with the mean flow should be Cr = Us ■ However, the wave speed 
Cr 7^ Us holds for — //i < ^^^^ < , as has been pointed out in proposition 1. In short, the waves 
have no interaction with the mean flows whatever there is an inflection point or not. This conclusion 
can be easily generated to the flows in (3 plane. Thus, the waves have no interaction with the stable 
flows. This is the reason why the disturbed waves can't take energy from the stably mean flow. 
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On the other hand, [15] pointed out that of a unstable wave must lie between the minimum 
and the maximum values of the mean velocity profile. Thus, the unstable wave is stationary relative 
to the mean flow, and can take energy from the mean flow. So there are wave-mean flow interactions 
in unstable flows. 

§ 1.1.4 Conclusion 

In summary, the general stability criteria are obtained for inviscid parallel flow. Both the criteria 
and the proofs are remarkably simple and easy to understand, comparing to Arnol'd's nonlinear 
theorems. The new criteria extend the former theorems proved by Rayleigh, Tollmien and Fj0rtoft. 
This may shed light on the flow control and investigation of the vortex dynamics. Based on the 
stability criteria, the reason why the disturbed waves can't take energy from the stably mean flow is 
explained that there is no wave-mean flow interaction at all. 

§1.2 Longwave Instability 

§ 1.2.1 Introduction 

Shear instability, caused by the velocity shear, is one of most important instabilities in the flows. 
Although the mechanism of shear instability is not full revealed yet, it is applied to explain the 
instability of mixing layer, jets in pipe, wakes behind cylinder, etc. Some simple models were employed 
to study the shear instability, including Kelvin-Helmholtz (K-H) model, piecewise linear velocity profile 
[16] and continued velocity profiles [1], etc. To reveal the mechanism of shear instability, the stability 
of perturbation waves should be understood. In K-H model, the growth rate of the disturbance wave 
is proportion to the product of wavenumber and velocity shear, thus the shortwaves are more unstable 
than longwaves in K-H model. However, Rayleigh [16] found that the piecewise linear profile is linearly 
unstable only in a finite range of wavenumbers < k < kc , which means the shortwaves are stable 
in this case [2, 7, 16, 17]. A contradiction [2, 7] emerges as shortwave and longwave dominate the 
shear instability in K-H model and piecewise linear profiles, respectively. This was explained either as 
viscous effect must be considered for shortwave [2] or as longwaves do not "feel" the finite thickness 
of the layer [2, 7]. 

To dispel the contradiction, we should find out which dominates the shear instability, shortwave 
or longwave? It is from the numerical simulation (see [2, 7, 13] and references therein) and some 
theoretical analysis (e.g. Tollmien [5] and Lin [18]) that longwave might be more unstable than 
shortwave. They also proved the long-wave instability subject to the velocity profile U is either 
symmetric or monotone. We will investigate this problem following the way by Sun [19]. 

Moreover, a prior estimation of growth rate, which may be obtained from the investigation, is 
useful for unstable flows. For example, Howard's semicircle theorem [15] has been used to validate 
the numerical calculations [13]. Another useful result was found by H0iland [20] and Howard [15] that 
the growth rate uJi must less equate to half of the maxim of vorticity, i.e., < |?7'|max/2 • But this 
estimation is too cursory to use in applications. For example, U' is always great than zero even the 
velocity profile has no inflection point. So this estimation is trivial for these cases. Here we will show 
a refinement estimation of growth rate, which can be applied to general velocity profiles in parallel 
flows. 
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The motivation of this short letter is to investigate these problems within the context of inviscid 
parallel flow. The aim here is to find out some general characters for unstable waves. 



§1.2.2 Model and Analysis 

To this purpose, long- wave instability in shear flows is investigated via Rayleigh's equation [1, 2, 
7, 17]. For a parallel flow with mean velocity U{y) , where y is the cross-stream coordinate. The 
streamfunction of the disturbance expands as series of waves (normal modes) with real wavenumber 
k and complex frequency to = LOr + ioJi , where LOi relates to the grow rate of the waves. The flow is 
unstable if and only if uji > . We study the stability of the disturbances by investigating the growth 
rate of the waves, this method is known as normal mode method. The amplitude of waves, namely 
4> , holds 

- k^ct>) - -^<A = 0, (1-13) 
U — c 

where c = u/k = Cr + ici is the complex phase speed. The real part of complex phase speed Cr = oJr/k 
is the wave phase speed. This equation is to be solved subject to homogeneous boundary conditions 

(P = Oaty = a,b. (1-14) 
From Rayleigh's equation, we get the following equations: 



b 



"' + k^n') + ^''^^_"2^ M']dy = 0, (1-15) 



and 



Before the further discussion, we need estimate the rate of /j* Ijc/i'lpdy to ||(/>|pciy , as Sun did 
[19]. This is known as Poincare's problem: 

\<P'fdy = fi fufdy, (1-17) 

J a 

where the eigenvalue fi is positive definition for (p ^ . The smallest eigenvalue value, namely , 
can be estimated as fii > (jr^)^ • And an auxiliary function f{y) = J^jj is also introduced, where 
f{y) is finite at inflection point. 

With the preparation above, we have such consequence. If —f{y) < Q < oo , where Q is a 
positive constant, then the disturbances with shortwaves k > k^ are always stable, where kc is a 
critical wavenumber subject to k1 = Q — . We will prove the consequence by two steps. At first, 
we prove proposition 1: if Cr = Ug , the disturbances with shortwaves k > k^ are always stable. 

Proof: Since U" = f{y){U — Us) and Cr = Us , this yields to 
\\U — c r U — U^ 
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and 

U — c r 



/ 

•J a 



b 



(1-19) 



This contradicts Eq.(l-15). So proposition 1 is proved. 

Then, we prove proposition 2: if Cr ^ Us , there must be c? = with k > kc- 

Proof: Otherwise if cf 7^ , so according to Eq.(l-15) and Eq.(l-16), for any arbitrary real 
number Ut which does not depend on y , it holds 



' k'm') + K^-^m']dy = o. (1-20) 



But the above Eq.(l-20) can not be hold for some special Ut ■ For example, let Ut = 2cr — Us , then 
there is {U - Us){U -Ut)<\\U - cf , and 

\\U — c\\'^ \\U — c||^ 

For k > kc , this yields to 



' U "{U - Ut ) 

U-cW^ 



m'f + [e + :;' ]\\mdy > o, (1-22) 



which also contradicts Eq.(l-20). So proposition 2 is also proved. These two propositions are natural 
generalization of stabile criterion proved by Sun [19]. 

From the above two propositions, we can draw a conclusion that the disturbances with shortwaves 
k > kc are always stable. This means that the shear instability in flows must be long-wave instability. 
Furthermore, the short-wave stability means that without any viscous effect, the shortwaves can also 
be damped by shear flow. This mechanism is unlike the viscous mechanism that the viscosity has 
a damping effect on especially the shortwaves. It implies that the shear flow itself can damp the 
shortwaves. This is very general and important conclusion, which explains why the instabilities found 
in shear flows are mostly long-wave instabilities. 



§ 1.2.3 Discussion 

This result is also very important for numerical calculation, which means shortwaves can be 
truncated in the calculations without changes the stability of shear flow. So the growth rates of 
longwaves have enough information for judging the stability of shear flow. On the other hand, the 
truncation of longwaves would probably change the instability of the shear flow. So the streamwise 
length scale must be longer enough to have longwaves for the numerical simulations in shear flows, 
such as plane parallel flow and pipe flow. Otherwise the instability of shear flow would be damped 
without long- wave perturbations. 

Then the growth rate of unstable waves can be estimated here by following the former investiga- 
tion. Similar to the assumption above, —f{y) < Q = p'^fJ-i < oo , where 1 < p < 2 , we will show that 
the growth rate is subject to u>i < {p — I) y/JnlUtnax — t^minl ■ 
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Proof: It is from Eq.(l-15) that gives 

^k'm'dy= ['-[UT + ^Q^^Uf]dy. (1-23) 
Substituting Eq.(l-17) to Eq.(l-23) and recalling that fj^i < fj, , this yields 

' k'm' dy<l' - [^^^^ + f^im' dy. (1-24) 
Multiplying above inequality (1-24) by c? , we get 

u^f fufdy< f h{y)Ufdy, (1-25) 

J a J a 

where 

%) = -[^9^^ + /^i]c?- (1-26) 

Suppose the maxim of h{y) is , then the growth rate is subject to 

ujf < (1-27) 



hiy) < -[ ,,^ ^:,,7;% + Mi]c?. (1-28) 



Hy) < fii r .r/\ ,2l:2 M cf. (1-29) 



This follows 

U"{U 
-(Z^)2 + cf 

Substitution of /(y) into Eq.(l-28), gives 
.[p' - l)(U - c,)^ - cl 2 

When 

= {p - 1){U - Cr)^ (1-30) 
the right hand of Eq.(l-29) get its largest value 

P^ = {p-lfl2l{U-Crf. (1-31) 

Then the growth rate must be subject to 

(^i<{p- 1)\/W|t^max - f^iinl, (1-32) 

where C/min and C/max are minimum and maximum of U (?/) , respectively. And the wavenumber A^^iax 
corresponding to the largest growth rate is 



A;max = ViP- l)/^i- (1-33) 

So the result is proved. One should note that the fast growth rate coi is only an approximation, but 
not a precision one, so as to wavenumber A;max • 

Comparing with the previous one by H0iland and Howard, the above estimation includes not 
only a more precision estimation about the the growth rate of unstable waves but also the regime of 
unstable wavenumbers. This would be much helpful for validation in numerical calculations. 

As well known, the instability due to velocity shear is always associated to Kelvin- Helmholtz 
instability (K-H instability), in which the disturbances of all wavelengths are amplified. According 
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to K-H instability, the shorter the wavelength is, the faster the perturbation wave amplifies. As the 
above investigation points out that the shortwaves are always more unstable than longwaves in this 
case, the contradiction at the beginning can be dispelled. 

An physical explanation [2, 7] to the contradiction is that the K-H instability model has no 
intrinsic length scale, while Rayleigh's model has width of shear layer as length scale of waves. This 
can be noted from that Rayleigh's case reduces to the Kelvein-Helmholtz vortex sheet model in the 
long- wave limit A; ^ 1 [2, 7], which is dramatically explained as longwave do not "feel" the finite 
thickness of the shear layer [7]. Here we will show that this explanation can be extend to shear 
flows. It is from Eq.(l-33) that the fastest growing wavenumber fcmax is proportion to ^/JIi , then the 
corresponding wave length Amax is approximately 2(b — a)l yj]) — 1 . So the thinner the shear layer 
is, the lager the fastest growing wavenumber is. It can be seen that the asymptotic case of infinite 
small shear layer leads to K-H instability. This is another evidence that K-H instability is in essence 
of long-wave instability. In this case, the K-H instability is an approximation of shear instability when 
the wave length of perturbation is much longer than the width of shear layer. 

In summary, two general properties of shear instability are obtained in the investigations. Firstly, 
the shortwaves are always stable in the continued profile flows and the shear instability is due to 
long-wave instability. This implies that the shear flow itself can damp the shortwaves and that K-H 
instability is in essence of long-wave instability in shear flow. The result can be used to understand 
the phenomena in hydrodynamics instabilities. Secondly, the growth rate of unstable flow can be 
estimated as tOi < {p — l)y^(t^max — Umm) ■ This estimation extend the previous result obtained by 
H0iland and Howard. Both results are important in numerical applications. The first one provides 
the estimation of unstable wavenumbers, and the second one provides the estimation of growth rate 
of unstable waves. These results may be useful on both numerical calculation and stability analysis. 



§1.3 Mechanism of Shear InstabiUty 

§ 1.3.1 Introduction 

The hydrodynamic instability is a fundamental problem in many fields, such as fluid dynamics, 
astrodynamics, oceanography, meteorology, etc. There are many kinds of hydrodynamic instabilities, 
e.g., shear instability due to velocity shear, thermal instability due to heating, viscous instability 
due to viscosity and centrifugal instability due to rotation, etc. Among them, the shear instability 
(inviscid) is the most important and the simplest one, which has been intensively explored (see [2, 6, 7] 
and references therein). Both linear and nonlinear stabilities of shear flow have been considered, and 
some important conclusions have also be obtained from the investigations. 

On the one hand, the nonlinear stability of shear flow has been investigated via variational 
principles. Kelvin [21] and Arnol'd [9, 10, 14] have developed variational principles for two-dimensional 
inviscid flow [22] . They showed that the steady flows are the stationary solutions of the energy H . 
And if the second variation 5^H is definite, then the steady flow is nonlinearly stable. Moreover, 
Arnol'd proved two nonlinear stability criteria, and that the flow is linearly stable as 6^H is positive 
definite [10, 14, 22, 23]. However, Arnol'd's theorem is not convenient to use as 6^H is always 
indefinite in sign, except for two special cases (see [23] and references therein). Due to the lack of the 
explicit expressions in both H and the stability criteria the variational principle is inconvenient for 
real applications. 
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On the other hand, the hnear stabihty of shear flow has also been investigated via Rayleigh's 
equation. Within the linear context, there are three important general stability criteria, namely 
Rayleigh's criterion [1], Fj0rtoft's criterion [4] and Sun's criterion [19]. As all the criteria have explicit 
expressions, they are more convenient in real applications, and are widely used in many fields. Based 
on the previous investigations, Sun [19] also pointed out that the flow is stable for Rayleigh's quotient 
/(/) > (see Eq.(l-ll) behind). However, the sufficient criterion for instability still lacks. 

To understand the shear instability, some mechanisms were suggested. Among them, Kelvin- 
Helmholtz instability (K-H instability) is always taken as a prototype, which is physically explained 
as the instability of a sheet vortex [24] . An another mechanism of instability is due to the resonance 
of waves [2, 25, 26, 27, 28]. Butler and Farrell [26] clearly showed with numerical simulations that the 
resonance introduces an algebraic growth term into the temporal development of a disturbance. Baines 
and Mitsudera [27] also used broken-line profile velocity as a prototype to explain the interaction of 
waves. Their explanation is so brilliant that it can explain why the instability occurs for a finite range 
of wavenumber and how the waves amplify each other. However, both mechanisms are independent 
of the basic flows, for that Kelvin-Helmholtz model deals only with vortex and resonance mechanism 
only considers the waves. Thus the relationships between those mechanisms and the shear flows are 
still covered. 

Overview the previous investigations, the essence of the shear instability is remained to be elu- 
cidated. Yet both the explicit expression and the meaning of H have not been revealed before, 
although it is very important in variational principles. The connection between linear and nonlinear 
stability criteria should be retrieved explicitly. The relationships between instability mechanisms and 
the physical explanation for shear instability are needed. The aim of this paper is to reveal the essence 
of the shear instability by investigating the inviscid shear flows in a channel. And this would lead to 
a more comprehensive understanding on shear instabilities. 



§ 1.3.2 Stable Criterion 

For the two-dimensional inviscid flows with the velocity of U , the vorticity ^ = Vx [7 is conserved 
along pathlines [7, 22, 24]: 

| = | + (t,.V« = 0. (1-34) 

Its linear disturbance reduces to Rayleigh's equation provided the basic flow U being parallel. Con- 
siders an shear flow with parallel horizontal velocity U{y) in a channel, as shown in Fig. 1-5. The 
amplitude of disturbed flow streamfunction -0 , namely (j) i satisfies [2, 6, 7] : 

(0" - feV) - = 0, (1-35) 

where k is the nonnegative real wavenumber and c = Cr + ici is the complex phase speed and double 
prime " denotes the second derivative with respect to y . The real part Cr is the phase speed of 
wave, and Q 7^ denotes instability. This equation is to be solved subject to homogeneous boundary 
conditions (p = at y = a,b . 

Considers that the velocity profile U{y) has an inflection point yg at which U'J = U"{ys) = 
and Us = U{ys) ■ As Sun [19] has pointed out, the following Rayleigh's quotient /(/) > implies the 
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flow is stable. 

Iif)=mm— 7 l-oD 

where f{y) = j/^jj^ ■ While if /(/) < , then there is a neutral stable model with kjf = —I{f) and 
Cr = Us . Moreover, there are unstable modes with Cr = Ug and Q 7^ if /(/) < and < k < ■ 
This can also be proved by following the way by Tollmien [5], Priedrichs [6, 29, 30] and Lin [18]. Thus 
/(/) = means the flow is neutrally stable. And there is only one neutral mode with k = and 
Cr = Us in the flow. These conclusions can be summarized as a new theorem. 

Tollmien-Fridrichs-Lin theorem: The flows are neutrally stable, if /(/) = . The flows are stable 
and unstable for /(/) > and /(/) < , respectively. 

From the above proof, it is obvious that the shortwaves (e.g. k ^ k^) are always more stable 
than longwaves (e.g. k <C k^ ) in the inviscid shear flows [30, 31]. So the shear instability is due to 
long-wave instability, and the disturbances of shortwaves can be damped by the shear itself without 
any viscosity [31]. 

As mentioned above, the linear stability criterion can be derived from the nonlinear one [14]. 
Moreover, the nonlinear stability criterion can also be obtained from the linear one. To illuminate 
this, the nonlinear criterion is retrieved explicitly via the above theorem, which is briefly proved as 
follows. 

Similar to Arnol'd's definition, the general energy H here is defined as 

1 /■\1„„.,,2 



^=2]^ i^W^W + hmdy (1-37) 

where ^(y) is the streamfunction of the fiow with ^ = U{y) — Us , and h{'^) is a function of ^ . 
The variation of 5H = gives 

A^- = h\^). (1-38) 

So Arnol'd's variational principle is retrieved. And the function h = ||V^|p/2 can also be obtained 
by solving Eq.(l-38). Then, the second variation d^H holds 

6^H = i ^'[11 VVf + h"m^]dy, (1-39) 

where ip denotes the variation 6^ . Noting that h"{^) remains unknown, Arnol'd's nonlinear criteria 
have not be extensively used. Fortunately, we can obtain the explicit expression of h"{^) here via 
the investigation on the linear stability criterion. For ^ = U — Us , 

^dJl^dh^d^^ U" ^ 

d^> dyd^ {U-Us) ^ ' 

As h"{^) is solved explicitly, 6^H has an explicit expression, which is greatly helpful for real appli- 
cations. 

Let the streamfunction of the perturbation expands as travelling waves il^{x, y, t) = cj){y)e^'^^^~^^'^ , 
where w is the frequency. The averages of ||V'i/'|P and HV'lP are (||(/)'|p + /c^||(/)|p)/2 and ||0|P/2 
along the fiow direction x , respectively. So Eq.(l-39) reduces to 

5^H = I [\\\<p'f + k^Uf + fiy)Uf]dy. (1-41) 
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The sign of 6^H is then associated with /(/) in Eq.(l-36). If /(/) < 0, the second variation 6^H 
can be both negative and positive, i.e., the stationary solution is a saddle point. And /(/) > 
implies 5^H is positive definite and vice versa. So the stable flow has the minimum value of the total 
kinetic energy H . The physical meaning of H can also be revealed, as the explicit expressions of H , 
^ = U-Us and h"{^) = f{y) have been obtained. 

First, according to the expressions, the velocity U in vorticity conservation law Eq.(l-34) can be 
decomposed to two parts: the rotational flow U — Us and the irrotational advection flow Us ■ The 
vorticity ^ inEq.(l-34) depends only on [/ — J7s,and Ug is only advection velocity. Then U — Ug and 
Us are associated with the dynamics and kinetics of the flow, respectively. Eq.(l-34) physically shows 
that the vorticity field is advected by Us , which can also known from the conservation of vorticity in 
the inviscid flows. A similar example is the dynamics of vortex in the wake behind cylinder, where the 
vortices dominate the dynamics of the flow and they are advected by mean flow (see Fig. 2 in [32]). 
The decomposition of velocity may be useful in vortex dynamics, for that our investigation clearly 
shows that the dynamics of the flow is dominated by vorticity distribution. 

Then the physical meaning of H can also be understood from the above investigation. It is not 
U but U — Us that is associated with the general energy H ,so H = \ £{U - Usfdy is not the 
total kinetic energy but the kinetic energy of flow with vorticity. Thus the stable steady states are 
always minimizing the kinetic energy of the flow associated with vorticity. This is also the reason why 
the How with maximum vorticity might be unstable, as Fj0rtoft's criterion shows. We would like to 
restate it as a theorem and to name it after Kelvin [21] and Arnol'd [14] for their contributions on 
this field [22]. 

Kelvin-Arnol'd theorem: the stable flow minimizes the kinetic energy of flow associated with 
vorticity. 

Both Tollmien-Fridrichs-Lin theorem and Kelvin-Arnol'd theorem are equivalent to the following 
simple principle [33]: The flow is stable, if and only if all the disturbances with Cr = Us are neutrally 
stable. 

§ 1.3.3 Physical explanation 

We have obtained the sufficient and necessary conditions for instability, then the physical mech- 
anism of instability can be understood from them. According to Fj0rtoft's and Sun's criteria [19], the 
necessary conditions for instability require that the base vorticity must be concentrated enough (e.g. 
sheet vortex). We call it "concentrated vortex" for latter convenience. The following investigation 
will reveal that the essence of shear instability is due to the interaction between the "concentrated 
vortex" and the corresponding resonant waves. 

As mentioned above, the concentrated vortex is a general model of sheet vortex in the Kelvin- 
Helmholtz model, for that the sheet vortex can be recovered as the concentrated vortex ^(i/s) oo (see 
[12] for a comprehensive discussion about the Kelvin-Helmholtz model and continued shear profiles). 
Then how the shear flow becomes unstable, if there is a concentrated vortex? As the sufficient 
condition for instability is /(/) < , the normal modes in the regime of < /c < ki\f with Cr = Us 
are unstable. These modes are stationary or standing waves, comparing to the velocity at inflection 
point. So that the shear instability is due to the disturbance of concentrated vortex by the standing 
waves with Cr = Us ■ In this case, the resonance mechanism is valid. The interaction waves propagate 
at the same speed with the concentrated vortex, so that they are locked together and amplified by 
simple advection [27]. In short, the disturbances on concentrated vortex is amplified like that in 
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Kelvin-Helmholtz model. 

This instability mechanism combines both K-H instability and resonance mechanism. Physically, 
the standing waves (with Cr = Us) can interact with the concentrated vortex, so they can trigger 
instability in the flows. While the travelling waves (with Cr ^ Us) have no interaction with the 
concentrated vortex, so that they can not trigger instability in the flows. This is the mechanism of 
shear instability. As pointed out above, the inviscid shear instability is due to long-wave instability. 
If the longwaves are unstable, they can obtain the energy from background flows. Thus, the energy 
within small scales transfer to and concentrate on large scales. In a word, the shear instability itself 
provides a mechanism to inverse energy cascade and to maintain the large structures or coherent 
structures in the complex flows. 

To illuminate the mechanism of shear instability, a physical diagram is also presented here by 
following the way of interpreting the K-H instability [24]. Fig. 1-5 sketches the mechanism of shear 
instability in terms of wave disturbances of vortices. The mean velocity profile U{y) has an inflection 
point at Us = with Us = , and the corresponding vorticity is ^(y) = —U'{y) . There is a local 
maximum at ys in the unstable vorticity profile in Fig.l-5a. According to Eq.(l-34), the vorticity is 
conserved in the inviscid flows. If the vortices at the local maximum (A, B and C) are sinusoidally 
disturbed from their original positions (dashed line) to new places (solid curve), they have negative 
vorticities with respect to the undisturbed ones. The vortices will induce cyclone flows around them 
in consequence. The flows around the vortices become faster (slower) in the upper (lower) of vortex A 
referred to the basic flow U (y) . The pressures at upper and lower decrease (indicated by - signs) and 
increase (indicated by + signs) according to Bernoulli's theorem, respectively. Then vortex A gets 
a upward acceleration due to the disturbed pressure difference as the uparrow shows. This tends to 
take the vortex away from its original position, so the flow is unstable. On the other hand, Fig.l-5b 
depicts the disturbances in a stable velocity proflle, where a local minimum is in the vorticity proflle. 
The disturbed vortices have positive vorticities with respect to the undisturbed ones. The vortices 




Figure 1-5: Sketch of shear instability: physical interpretation. Left parts depict the proflles of velocity 
U{y) and vorticity ^(y) , right ones depict the disturbance of vorticities. The unstable veloctiy profile 
U{y) has a local maximum in vorticity ^(y) (a). If the vortices (A, B and C) disturbed from their 
original positions (dashed line) to new places (solid curve), they will be taken away from their original 
positions due to pressure difference. The stable veloctiy profile U {y) has a local minimum in vorticity 
i{y) (b)- The disturbed vortices (A', B' and C) will be brought back to their original positions due 
to pressure difference. 
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will induce anticyclone flows around them in consequence. So vortex A' get a downward acceleration 
due to the disturbed pressure difference as the downarrow shows. This tends to bring the vortex 
back from its original position, so the flow is stable. In this interpretation, the advection of Ug is 
independent of the shear instability, only the fiow field of U — Us and the corresponding vorticity ^ 
are the dominations. The unstable disturbances in Fig. 1-5 have Cr = Us , which consists with the 
above discussions. This physically explains why the maximum and minimum vorticities have different 
stable aspects. 

Though Tollmien-Pridrichs-Lin theorem is a sufficient and necessary condition for stability, the 
unknown in the theorem restricts its application. So some simple criteria may be more useful. For 
the parallel flows within interval a < y < b , there are two simple criteria. 

Corollary 1: The flow is stable for f{y) > -{^f [19, 33]. 

Corollary 2: The flow is unstable for f{y) <-{^f [33]. 

In summary, the general stability and instability criteria are retrieved for inviscid parallel flow 
within linear context, which are associated with the nonlinear stability criteria, i.e., minimizing the 
kinetic energy of flow. Then the mechanism of shear instability is explained as the resonance of 
standing waves with the concentrated vortex at Cr = Ug ■ The physical process is also sketched by 
extending the way of interpreting K-H instability. Finally, some useful criteria are given. These results 
would lead to a more comprehensive understanding on shear instabilities, especially for undergraduate 
students. 

§1.4 Rotating Flow 

§ 1.4.1 Introduction 

The instability of the rotating flows is one of the most attractive problems in many fields, such 
as fluid dynamics, astrophysical hydrodynamics, oceanography, meteorology, etc. Among them, the 
simplest one is the instability of pure rotation flow between coaxial cylinders, i.e., Rayleigh- Taylor 
problem, which has intensively been explored [6, 17]. 

Two kinds of instabilities in inviscid rotating flow have been theoretically addressed in the liter- 
atures. One is centrifugal instability, which was first investigated by Rayleigh [1, 6]. He derived the 
circulation criterion for the inviscid rotating flows that a necessary and sufficient condition for stability 
to axisymmetric disturbances is that the square of the circulation does not decrease anywhere. This 
criterion is also be stated as the Rayleigh discriminant $ > (see Eq.(l-46) behind), and is always be 
used in astrophysical hydrodynamics [34]. It is also generalized to non-axisymmetric fiows [35]. The 
other is known as instability due to two-dimensional disturbances in rotating flows, which is similar 
to the shear instability in parallel flow. We call this instability as the shear instability in rotating 
flow hereafter. For this instability, Rayleigh also obtained a criterion, i.e., inflection point theorem in 
inviscid rotating flows, which is the analogue of the theorem in parallel flows [1]. Following this way, 
Howard and Gupta [36] found a stability criterion for two-dimensional disturbance in inviscid rotating 
flow. However, the theoretical results remains scarce, due to the complexity of rotating flow. 

Comparing this to shear instability in parallel flows, the criteria for parallel flows are much more 
abundant. Fj0rtoft [4] and Sun [19] proved some more strict criteria. And the stability of two- 
dimensional in a rotating frame was also addressed by Pedley [37], which seems to be more complex 
than the stability problem in the pure rotation flows. 
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Motivated, then, by the theoretical criteria for parahel flows [4, 19], our study focuses on the 
instability due to shear in inviscid rotating flows. The aim of this letter is to obtain such criteria for 
the inviscid rotating flows, and the relationship between previous criteria is also discussed. Thus other 
instabilities may be understood via the investigation here. 

§ 1.4.2 Howard-Gupta Equation 

For this purpose, Howard-Gupta equation (hereafter H-G equation) [36] is employed. To obtain 
H-G equation, Euler's equations [2, 6, 17, 24] for incompressible barotropic flow in cylindrical polar 
coordinates (r, 6) are then given by 

9«r _^ ^ 9wr _^ _ _ _1 9p (142) 

dt dr r do r p dr^ 

and 

du0 du0 uedue , n^ug 1 dp 

^ + nr.— + — H = — . 1-43 

ot or r ou r pr od 

Under the condition of incompressible barotropic flow, the evolution equation for the vorticity can be 
obtained from Eq.(l-42) and Eq.(l-43), 

where ^ = ^-§p{rue) — is the vorticity of the background flow. Eq.(l-44) can also be derived 

from Fridman's vortex dynamics equation [22, 24]. And it admits a steady basic solution, 

= 0,ug = V{r) = 0(r)r, (1-45) 

where il(r) is the mean angular velocity. And Rayleigh discriminant is defined by 

$ = L±(nry. (1-46) 

§ 1.4.3 Stable Criterion 

Then, consider the evolution of two-dimensional disturbances. The disturbances ip'{r, 6, t) , which 
depend only on r , ^ and t , expand as series of waves, 

^'{r,e,t) = (l){r)e'^'''^-'^'\ (1-47) 

where (j){r) is the amplitude of disturbance, n is real wavenumber and uj = uJr + iooi is complex 
frequency. Unlike the wavenumber in Rayleigh's equation for inviscid parallel flows, the wavenumber 
n here must be integer for the periodic condition of 9 . The flow is unstable if and only if cjj > . In 
this way, the amplitude (j) satisfies 

2 

(nn-uj)\D*D-^](t>--{Di)4> = Q, (1-48) 

where D = d/dr , D* = d/dr + l/r . This equation is known as H-G equation and to be solved subject 
to homogeneous boundary conditions 

D(j) = at r = ri,r2. (1-49) 
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By multiplying J_^^^ to H-G equation Eq.(l-48), where cp* is the complex conjugate of cp , and 
integrating over the domain ri < r < r2 , we get the following equation 

r{cP%D*D)4> - + = 0. (1-50) 

Then the integration gives 

r{\Wf + [4 + "?"""?f ]ll^f}d^ = 0, (1-51) 



where (j)' = D(j) , ^' = -D(^) and iv* is the complex conjugate of lo . Thus the real part and image 
part are 



- (0 - Cr)^' 



-iUT + G+ '^^:7T^m'}dr = 0, (1-52) 



and 



"-"2^r=0, (1-53) 



117 - c 



12 



where c = to/n = Cr + ici is the complex angular phase speed. Rayleigh used only Eq.(l-53) to prove 
his theorem: The necessary condition for instability is that the gradient of the basic vorticity ^' must 
change sign at least once in the interval ri < r < r2 ■ The point at r = is called the inflection 
point with = , at which the angular velocity of = ^{fs) ■ This theorem is the analogue of 
Rayleigh's inflection point theorem for parallel flow [1, 6]. 

Similar to the proof of Fj0rtoft theorem [4] in the parallel flow, we can prove the following 
criterion. 

Theorem 1: A necessary condition for instability is that — ^s) < (or — 17^) < 0) 

somewhere in the flow field. 

The proof of Theorem 1 is trivial, and is omitted here. This criterion is more restrictive than 
Rayleigh's. Moreover, some more restrictive criteria may also be found, if we follow the way given by 
Sun [19]. If the velocity profile 0,{r) is stable (cj = 0), then the hypothesis Q 7^ should result in 
contradictions in some cases. So that a more restrictive criterion can be obtained. 

To find the criterion, we need estimate the rate of f^^ r||(/>'|pdr to f^^ \\(/)\\'^dr , 

\\\())'fdr = /u / ' Ufdr, (1-54) 

ri J ri 

where the eigenvalue is positive definition for (j) ^ . According to boundary condition Eq.(l-49), 
(j) can expand as Fourier series. So the smallest eigenvalue value, namely /Ui , can be estimated as 
m > ri7r2/(r2 - . 

Then there is a criterion for stability using relation (1-54), a new stability criterion may be found: 
the flow is stable if —(^1 + l/r2) < qSj^ < everywhere. 

To get this criterion, we introduce an auxiliary function /(r) = ^ , where /(r) is finite at 
inflection point. We will prove the criterion by two steps. At first, we prove proposition 1: if the 
velocity profile is subject to — (/Ui + l/r2) < f{r) < , then 7^ l^s • 
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Proof: Since — (/ii + l/?"2) < fi'i') < , then 
and if Cr = 0,s and 1 < n , this yields to 



(1-56) 



This contradicts Eq.(l-52). So proposition 1 is proved. 

Then, we prove proposition 2: if — (/ii + l/r2) < /(r) < and c,. 7^ fig , there must be cf = . 

Proof: If c? 7^ , then multiplying Eq.(l-53) by (c^ — Ct)/ci , where the arbitrary real number q 
does not depend on r , and adding the result to Eq.(l-52), it satisfies 



r\\Q — c"^ 



riur + B + ::i m'} dr = o. (1-57) 

n 



But the above Eq.(l-57) can not be hold for some special q . For example, let q = 2cr — 0,s , then 
there is {Q, — 0,s){Q — q) < ||Q — c|p , and 

This yields 

'[rll'A'f + (- + ^-^^^)\m']dr > 0, (1-59) 

ri II' ' ~ C|p 

which also contradicts Eq.(l-57). So the second proposition is also proved. 

Using 'proposition 1: if — (/^i + 1/^2) < f{r) < then Cr / ' and 'proposition 2: if —{fii + 
1/^2) < /(r) < and Cr / then q = ', we find a stability criterion. 

Theorem 2: If the velocity profile satisfy — (/xi + l/'r2) < /(r) < everywhere in the flow, it is 
stable. 

This criterion is the analogue of the theorem proved by Sun [19]. Both theorem 1 and theorem 
2 here are more restrictive than Rayleigh's theorem for the inviscid rotating flows. The theorems 
indicate the probability that a vorticity profile with local maximum or minimum would be stable, if 
it satisfies the stable criteria. Theorem 2 implies that the rotating flow is stable, if the distribution 
of vorticity is relatively smooth. As shown by Sun [19], the instability of inviscid parallel flows must 
have vortices concentrated enough. This is also the shear instability in rotating flows. Since several 
stable criteria for inviscid rotating flows have been obtained, it is convenient to explore the relationship 
among them, as discussed followed. 

The criteria for rotating flow can be applied to parallel flows, given narrow-gap approximation. 
First, Pedley's criterion is covered by the centrifugal instability criteria. As mentioned above, Pedley 
[37] considered the stability of two-dimensional flows [/ in a frame rotating with angular velocity . A 
criterion is found that instability occurs locally when 2^}{2Q — U') < , where U' = dU/dr represents 
radial shear of horizontal velocity. Pedley's criterion, which is recovered by later researches [38, 



is in essence the special case of Rayleigh's circulation criterion, i.e., ^(il^r"^) < for instability. 
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Here the proof is briefly given. Considering the narrow-gap approximation r2 — ri = d <^ ri and 
the large radii approximation 1/ri — > in Rayleigh's circulation criterion, there is fi'r ~ — [/'. So 
$ = d{^l'^r^)/dr/r^ = 2J7(2J7 — U') < 0, which is exactly Pedley's criterion. Thus Pedley's criterion 
is covered by Rayleigh's circulation criterion. Second, the stable criteria for parallel flows, such as 
Rayleigh's theorem [1] and Fj0rtfot's theorem [4], can be derived from those for rotating flows, given 
the narrow-gap approximation r2 — ri = d <C ri and the large radii approximation 1/ri . 
Following this way, the results of the Taylor-Couette system can also be applied to the plane Couette 
system [40]. The proof is omitted here, as the approach is trivial too. As pointed out by Sun [19], all 
of the shear instability criteria for two-dimensional flows are the special cases of Arnol'd's nonlinear 
criteria [14], which are much more complex yet not widely used. In general, all the known stability 
criteria for parallel flows (even in a rotating frame) can be derived from the stability criteria for 
rotating flows. 



§ 1.4.4 Discussion 

As well known, these two kinds of instabilities are independent with each other. Howard [6, 36, 41] 
has given an example which is stable to axisymmetric disturbances but unstable to two-dimensional 
(shear) disturbances. Here another example is given that the flow is unstable to axisymmetric distur- 
bances but stable to two-dimensional disturbances according the stability criteria above. To illuminate 
this, a simple example is given as il(r) = 1 — r with the vorticity ^ = 2 — 3r . As shown in Fig. 1-6, 
both (solid line) and (dashed line) are plotted within the interval ri < r < r2 ■ It is from 
Rayleigh's inflection point criterion that the flow is always stable to two dimensional disturbances. 
While the flow is unstable for the axisymmetric disturbances, if2/3<ri<l or 2/3<r2<l. 
Then, a shear stable flow may be unstable due to centrifugal instability. Considering that there are 
three-dimensional disturbances, the flow might be unstable, even it is centrifugal stable and shear 
stable. This might be the reason why there is no general criterion known for H-G equation when 
non-axisymmetric disturbances are considered [36]. 

And the relationship between rotation and stratification is also noted. As "Rayleigh observed 
that there is an analogy between the stability of rotating flows and the stability of a stratified fluid 
at rest in a gravitational field" [6], the criteria between them can also be analogy with each other. 
For the stratified flow with vertical velocity shear, Howard [15] obtained a stability criterion, i.e., the 
flow is stable if N"^ > U'"^ /A , where N is known as Brunt- Vaisala frequency, U is horizontal velocity 

2 

1 
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-3 

"*0 0.5 1 1.5 2 

Figure 1-6: Angular velocity (solid line) and corresponding vorticity ^ (dashed line) versus r . 
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and U' = dU /dz represents vertical shear of horizontal velocity. This criterion is then an analogue 
of the following criterion by Howard and Gupta [36]: The rotating flow is stable for axisymmetric 
disturbances if Eq.(l-60) holds, 



where W is the axis velocity, and it vanishes in the pure rotating flows, as in Rayleigh's criteria. 
Comparing the two stability criteria, the analogues of $ and dW/dr are N"^ and V , respec- 
tively. Physically the rotating establishes a potential distribution along the ratio, which is denoted by 
Rayleigh discriminant <I> . Similarly, the stratification establishes a potential distribution along the 
gravity, which is denoted by Brunt- Vaisala frequency. And the analogue of the radial shear in rotating 
flow is right the vertical shear in stratified fiow. So the criteria between rotation and stratification 
can also be analogy with each other in this way. 

In short, the general stability criterion is obtained for inviscid rotating fiow. Then Pedley's 
cirterion is proved to be an special case of Rayleigh's criterion. These results extend Rayleigh's 
infiection point theorem for curved and rotating fiows, and they are analogues of the theorems proved 
by Fj0rtoft and Sun for the two-dimensional inviscid parallel fiows. 

§1.5 General Barotropic Wave 

Barotropic waves are widely existed in the fiow with shear. However, the general theory for 
barotropic waves still lacks though the stability of linear waves (e.g. Rossby wave) are always the 
most concerned. The theory of wave-mean fiow interaction points out that the fiow is unstable, the 
disturbed waves can grow by taking energy from the mean fiow. However, why the disturbance waves 
can't take energy from the mean fiow in the stable fiow is not known. Nevertheless, the general criterion 
can be great helpful to understand these problems. As Arnol'd's criteria deal with the total energy, 
they can not provide such kind of information. In this way, a simply linear criterion corresponding 
to Arnol'd's second nonlinear stability theorem is necessary for considering the stability of waves. It 
also provides a new way to investigate shear instability. The aim of this short paper is to find such a 
stability criterion from Rayleigh's equation in normal mode way, through which the shear instability 
may be understood. 

The above criteria would be helpful for understanding the wave-mean fiow interaction, especially 
the Rossby wave- mean fiow interaction in barotropic fiows. According to the stable criteria, the 
necessary condition for wave-mean fiow interaction can be obtained. And why the disturbed waves 
can't take energy from the mean fiow in the stable fiow is revealed. If the fiow is stable, there is no 
wave-mean fiow interaction at all. 

First, when the velocity profile has no infiection point, then the speeds of barotropic waves 
Cr — U ^ Q . According to Eq.(l-3), U — Cr > Q holds for U" < 0, i.e., the barotropic waves are 
always west-propagation relative to the mean fiow U if the vortex gradient is positive. This extends 
the west-propagation theory of Rossby waves. And C/ — < holds for U" > , i.e., the barotropic 
waves are always east-propagation relative to the mean fiow U if the vortex gradient is negative. 
Second, when the velocity profile has an infiection point U'^ = , the speed of most favorite wave 
which might have interaction with the mean fiow should be = Us ■ However, the wave speed 
Cr 7^ Us holds for — /ii < jjzijj- < , as has been pointed out in proposition 1. In short, the waves 
have no interaction with the mean fiows whatever there is an infiection point or not. This conclusion 




(1-60) 
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can be easily generated to the flows in [3 plane. Thus, the waves have no interaction with the stable 
flows. This is the reason why the disturbed waves can't take energy from the stably mean flow. 

On the other hand, [15] pointed out that c,. of a unstable wave must lie between the minimum 
and the maximum values of the mean velocity profile. Thus, the unstable wave is stationary relative 
to the mean flow, and can take energy from the mean flow. So there are wave-mean flow interactions 
in unstable flows. 

In summary, the general stability criterion is obtained for inviscid both parallel and rotating 
flows. Then Pedley's cirterion is proved to be an special case of Rayleigh's criterion. These results 
extend Rayleigh's inflection point theorem for curved and rotating flows, and they are analogues of the 
theorems proved by Rayleigh, Tollmien and Fj0rtoft for the two-dimensional inviscid parallel flows. 
Besides, this would intrigue future research on the mechanism of hydrodynamic instability. 
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§2.1 Establishment of Horizontal Convection 

§ 2.1.1 Introduction 

The abyssal ocean circulation is thought of an important energy conveyor belt, which has great 
impact to climate change [42, 43, 44]. As former investigators noted the fact that density water 
sinks by surface cooling at North Atlantic Ocean, this is so called density-driven flow. Horizontal 
convection, in which the flow is uneven heated at the horizontal surface, was taken as a model of such 
circulation. However, a novel idea emerges when the ocean energy balance is considered [45, 46, 47]. 
The horizontal convection become a key model to exam the theories. Two classes of flows are often 
used. The first one, named as two-basin forcing flows hereafter, has symmetric surface forcing like 
that of South and North Atlantic Ocean basins. The multiple equilibria and bifurcation phenomena 
of such flows are often discussed, especially in numerical way [48, 49, 50, 51]. The second one, named 
as one-basin forcing flows hereafter, has monotone forcing from one side to another like that of North 
Atlantic Ocean basin. This kind of flow is more conveniently used in the experiments [52, 53, 54]. 

Some general properties of horizontal convection have been obtained by theoretic studies. Above 
all, [52] found the 1/5-power laws of Rayleigh number Ra for flow strength ^'max and the Nusselt 
number Nu , which are consistent with experiments and numerical simulations. This 1/5-power law 
of Ra is generally valid no mater what the flow is steady or non-steady [e.g. 48, 52, 53, 54, 55, 56]. 
Thus it is useful to consider what the horizontal convection would be as Ra — > oo . [51] focused on the 
energy dissipation under the conditions of the viscosity and thermal diffusivity k being lowered 
to zero with Prandtl number Pr fixed. They proved that the horizontal convection is non-turbulent 
under certain definition of turbulence. Motivated by above results, [56] claimed Nu is bounded by 
Ra^^ as Ra^oo. 

On the other hand, numerical simulations and experiments were also used to find details of the 
flows. [52, 55] found by both experiment and numerical simulations that the flows are steady and 
stable in one-basin circulation. Comparing to Rossby's one-basin forcing, [48] studied the multiple 
equilibria under symmetric two-basin forcing numerically. They also claimed that the flow is always 
symmetric if only theromal or salt forcing at least for their parameter regimes. If both kinds of forcing 
are considered, there will be symmetry breaking and even unsteady flows [e.g. 48, 49, 50, 57]. For 
one-basin circulation, it was thought that horizontal convection must be steady and stable even for 
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Ra oo , according to Sandstorm's theorem [e.g. 46, 58, 59]. 

But this scenario seems to be violated in recent investigations. [51] studied the horizontal convec- 
tion under the symmetric two-basin forcing. The smaller the Prandtl number is, the more unstable 
the flow is. As the Prandtl number increases to 10, the circulation tends to be a shallow cell. They 
also obtained the thresholds for the transition from steady flows to unstable and steady flows. These 
unsteady flows then also ware reported in experiments [53]. However, such unsteady convective mix- 
ing and turbulent interior motion in the exmpriments are nonturbulent according to the definition of 
turbulence [51]. As the plume rises from bottom to top and through the full depth of the tank, this 
kind of fiows are referred as full-penetrating fiows. While [54] showed in their experiments a totally 
different result. The motion of the circulation, though being visible to the naked eye, is vanishingly 
small. The convection cells, shallow and near the heating surface, are steady and stable. This kind of 
flows are quite different from those in [53], and they were referred as partial-penetrating flows [54]. 

Though lots of numerical simulations on the horizontal convection, none of them have obtained 
such partial-penetrating flow and the onset of such flow needs to reveal yet. The main purpose of this 
paper is to investigate the partial-penetrating flow by numerical simulation, thus resulting in a more 
comprehensive view on this issue. The model and the scheme are sketched in §§2.1.2, in which a 
benchmark solution is included. The establishment of the circulation and the partial-penetrating cell 
of horizontal convection are depicted in § § 2.1.3, with the onset of the partial-penetrating fiow being 
discussed. 



§2.1.2 Model and Scheme 

We consider the the horizontal convection flows within the two-dimensional domain, and Boussi- 
nesq approximation is assumed to these flows. The horizontal (y) and vertical (z) regimes are 
< y < L and < z < D , respectively. Similar to [51], the depth D is taken as reference length 
scale and A = D/L denotes the aspect ratio. We use ^ = 1 in present work, which is consistent 
with the experiments by [54]. Taking account of nondivergence of velocity flled in Boussinesq approx- 
imation, the lagrangian streamfunction ^ and the corresponding vorticity to are introduced. The 
velocity 1? = {v,w) , where horizontal velocity v = ^ and vertical velocity w = , respectively. 
The governing equations [48, 49, 50, 51, 56] in vorticity-streamfunction formulation are 



^ + J(^',cu) =- Pr(V2c^ + Ra^) (2-lb) 
ot oy 

V^^r = -co (2-lc) 

where J(^',0) = — denotes the nonlinear advection term. There are two important di- 

mensionless parameter in Eq.(2-1), i.e. Rayleigh number Ra = ax^TgL^ /{ku) and Prandtl number 
Pr = f/fv, where g , ax, AT, L, k and i' are gravity acceleration, thermal expansion coefficient, 
surface temperature difference, length of horizontal domain, thermal diffusivity and kinematic viscos- 
ity, respectively. The surface buoyancy forcing is T = sin(|y) , so that only one-basin flows instead 
of symmetric two-basin flows can be obtained. 

There are two important quantity describing the circulation, i.e. the non-dimensional stream- 
function maximum and the non-dimensional heat flux. The non-dimensional streamfunction maxi- 
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mum 'I'max = ^max/ ^ ) where \I'max is maximum of the dimensional streamfunction. For the 
non-dimensional heat flux is defined as /t = dT/dz at the heated surface. 

The above Eq.(2-1) is solved with Arakawa scheme [e.g. 60, 61] and non-uniform grids. Comparing 
to the other schemes, Arakawa scheme is more accuracy but more expensive, and it has also been 
applied to horizontal convection flows at high Rayleigh number [62]. 

To validate the scheme, we calculate the nature convection problem with the resolution of 80 x 80 
meshes. Table 2-1 compares the present solutions with bench mark solutions from [63], who solved 
the problems using a pseudo-spectral Chebyshev algorithm and [64], who used fourth-order compact 
finite difference scheme. It is obvious that the our scheme is very accuracy for such flows. 

§2.1.3 Results 

1. the establishment of circulation 

Similar to the experiments by [54], we use A = 1 for the numeric simulations, the circulations 
are obtained for Ra > 10^. The dimensionless time t = t*/T , where t* and r = D'^ j k are the 
dimensional and the unit scaling times, respectively. A typical value of t = 1 is about 80 hours in the 
dimensional time, given D = 20 cm and k = 1.4 x 10~^ err? j s , which are approximate to the values 
used in the experiments by [53, 54]. According to the flow pattern, the establishment of the circulation 
can be divided into three stages: (1) startup of circulation, (2) damp of secondary circulation, (3) 
ampliflcation of primary circulation. The former two stages are relatively short, while the third one 
is very long. At the end of last stage, the circulation is fully established. 

To illuminate this, the circulation of Pr = 8 and Ka = 5 x 10^ is taken as an example. During 
the flrst stage, the circulation is established as soon as the surface forcing is superposed (Fig.2-la). 
The gradient of horizontal buoyancy drives the water like a lid, so that a very energetic primary 
circulation ( ^I'max > 56 ) generates (Fig.2-16), which was also observed by the experimentalist [54]. 
But this process is much more faster, and soon the secondary circulation emerges below the primary 
circulation at i = 3.2 x 10"'^ (Fig.2-lc). The secondary circulation becomes stronger and stronger, 
and the primary circulation becomes weaker and weaker. At the end of the stage, the primary 



author 


Ra 




*max 


'Umax(2/) 


^max(2^) 


Nu 


Present 




16.430 


16.863 


64.47(0.85) 


219.17 (0.038) 


8.828 


[63] 


10^ 


16.386 


16.811 


64.83(0.85) 


220.56 (0.038) 


8.822 


[64] 


10^ 


16.386 


16.811 


64.83(0.85) 


220.57 (0.038) 


8.825 


Present 


10^ 


29.586 


30.426 


146.05(0.886) 


687.17 (0.022) 


16.535 


[63] 


10^ 


29.361 


30.165 


148.59(0.879) 


699.17 (0.021) 


16.523 


[64] 


10^ 


29.356 


30.155 


148.57(0.879) 


699.17 (0.021) 


16.511 



Table 2-1: Comparison of the bench mark solutions from Le Quere and [64], Pr = 0.71 . '^mid ) ^max 
are the values in the midpoint and the maximum of streamfunction, respectively. Umax and f max are 
the maximum of horizontal and vertical velocity component on the vertical midplane x = 0.5 and 
y = 0.5 , together with its position y and x , respectively. And ISlu is average Nusselt number at the 
heated wall. The resolution is 80 x 80 meshes for present results. 
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Figure 2-1: Snapshots of the flow fields (streamfunction ^ ), with solid counter curves for ^' > and 
dashed counter lines for ^ < . 
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Figure 2-2: Snapshots of the temperature fields corresponding to these in Fig. 2-1, with counter inter- 
vals being 0.1. 
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circulation becomes partial-penetrating, the secondary circulation becomes full-penetrating (Fig.2- 
1/). Meanwhile, the heat conducts from top to bottom along the side walls (Fig. 2-2). 

Comparing to the first stage, the second one is longer, during which the secondary circulation 
damps and breaks into several weaker circulations (Fig. 2-3). In consequence, the smaller circula- 
tions emerges from the bottom one after another. And the flow field fulfills such partial-penetrating 
circulations, as observed by [54]. The damp of the secondary circulation is a puzzle for the experi- 
mentalists, as the corresponding temperature field was not well measured then. Here the temperature 
fields (Fig. 2-4) explains the reason why such full-penetrating circulation can't be maintained. The 
lack of driven forcing and the viscous friction are the reasons. As in Fig. 2-4, there is vanishingly small 
buoyancy gradient in the secondary circulation, so that it damps down due to viscosity and boundary 
friction, especially the circulation near the bottom. Meanwhile, the primary circulation covers the 
convection region, where the horizontal buoyancy gradient is remarkably large. 

The last stage is a long and slowly process of approaching to quasi-equilibrium state. The primary 
circulation is amplified in this stage, and the secondary circulations disappear from the bottom, like 
their emerging process (Fig. 2-6). It is notable that the primary circulation (e.g. the shadowed shallow 
circulation cell in Fig. 2-6) seldom changes during this stage, which implies that the flow near surface 
approaches to quasi-equilibrium state relatively faster. While the process to quasi-equilibrium state 
is very slowly near the bottom. Comparing the flow field and temperature field, it is clear that the 
flow field can't be steady until thermal conduction is balanced. 

The above stages have different time scales. The startup stage is the most fast stage, during 
with the circulations are established within t = 10^^ (several minutes in laboratorial time). The 





y V y 

Figure 2-3: Snapshots of the flow fields (streamfunction \I'), with solid counter curves for ^' > 
and dashed counter lines for ^' < . The counter intervals are 8 and -1 for ^' > and ^I* < , 
respectively. 
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Figure 2-4: Snapshots of the temperature fields corresponding to these in Fig. 2-3, with counter inter- 
vals being 0.1. 

second stage and the last state are time scales of t = 10~^ (several hours) and t = 2.5 (several days), 
respectively. 

In the following simulations, we use the total kinetic energy Ef^. and ^'max as indexes to oversee 
the evolution of the circulations, where Ek = §{v'^ + w^) /2dy dz . It is notable that the time ap- 
proaching to quasi-equilibrium state is very long (about one week) in our numerical simulations, but 
is is remarkably shorter in the experiments [53, 54]. According to the parameters used by [54], not all 
of the experiments take enough time, so that the circulations are not fully established in some cases. 
This is more serious in the experiments by [53], where the experiment time of 30 hours is much less 
than of 200 hours. 

The establishment time te of the circulation is defined as the time when d^/dt < O.OOl^'max 
holds for whole field. It is found that the partial-penetrating circulation is established very fast, but 
it takes a very long time for the flow to approach to equilibrium state. Fig.2-7 displays the time 
evolutions of total kinetic energy E^ and ^'max with time t . The circulation is established very fast, 
and both Ej. and ^max reach their 95% of the equilibrium values within t = 0.05 (4 hours). In 
addition. Fig. 2-76 agrees well with the experiments (e.g. Fig. 2 in [54]), and this fast establishment of 
the circulation was also noted by [53]. However the total time to get the quasi-equilibrium state is 
relatively longer, it takes at least = 2.5 (200 hours) for the flow being steady in present simulation. 

Moreover, the numerical simulation shows that tg = 2.5 is widely valid, at least for all the cases 
within the regime of 1 < Pr < 10 and 10^ < Ra < 10^^ . It also implies that the establishment of 
the circulation is governed by conduction process along the side walls, as te is only determined by 
D and k . The conclusion of = 2.5 is also potentially useful for both numerical simulations and 
laboratorial experiments. 
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Figure 2-6: Snapshots of the temperature fields corresponding to these in Fig. 2-5. 
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Figure 2-7: The time evolution of the flow field for total kinetic energy Ek (a) and maximum flow 
^'max (b) at Pr = 8 and Ra = 5 x 10^ . 

The fast established but slowly steadied flow is due to the horizontal evolution of buoyancy (or 
temperature). It is well known that the horizontal gradient of buoyancy drives the circulation, which 
can also be known from Eq.(2-16). Thus the fast establishment of the temperature gradient near the 
surface makes establishment of the circulation very fast. To illuminate this, the temperature fields 
of the circulation are shown in Fig. 2-6. The horizontal temperature gradient is relatively large near 
the top surface, but it is relatively to small to be vanished near the bottom. Comparing Fig.2-5a 
with Fig. 2-6 a, the main circulation is right within the zone where temperature gradient is remarkably 
large, so as to the figures at other times. As the time goes on, the main circulation becomes deeper 
and deeper with the downward propagation of buoyancy gradient. So the slowly heat conduction from 
the top surface to the bottom makes the establishment of the circulation near bottom to be a long 
time process. Hence, the depth D of the tank and the thermal diffusivity k both determinate the 
establishment time tg of the circulation. And it is from Fig. 2-7 that is 2.5 for the circulation of 
Pr = 8 and Ra = 5 x 10^ . 



2. power laws 

As the parameters in different ocean global circulation models (OGCM) varies in wide regime 
[e.g. 47, 65], it is very useful to know how the sensitivity of the results to the parameters. This is also 
useful when the experimental results are extrapolated to real ocean circulations [53] . 

Then, the power laws at different Prandtl numbers (e.g. Pr = 1 , Pr = 4 , Pr = 8 and Pr = 10 ) 
are calculated to investigate the sensitivities of Nu and ^'max to Pr . Fig.2-8a and Fig. 2-86 show 
the 1/5-power laws of ^'max and Nu , respectively. Noting that ^'max instead of ^'max and Rap = 
RaPr^/"^ instead of Ra are in Fig.2-8a, it implies that the flow is dominated by thermal diffusivity 
K, . The larger Prandtl number is, the stronger the flow is. All the direct numeric simulation (DNS) 

1/5 

data at Pr = 1 , Pr = 4 , Pr = 8 and Pr = 10 lie in the line of this power law, \I'max = 0.75Rap , 
where Rap = RaFr^^'^ . On the other hand, [54] fitted their experiments data into two different 1/5- 
power laws. It is amazing that more than half of the experimental data (14/25) by [54], of which 
the reduced acceleration due to gravity g' > 2.5 cm/s^ , lie in this line (right triangles in Fig.2-8a). 
Taking account that the aspect ratio in DNS is not exact the same as in their experiments, these data 
are consistent well with each others. Meanwhile, the Nusselt number Nu in DNS lie in two slightly 
displaced parallel lines, one for Pr = 1 and the other for Pr = 4 , Pr = 8 and Pr = 10 . It seems 
that ^'max is more sensitive to Prandtl number than Nusselt number does. 
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It is notable that 1/5-power law in Fig.2-86 is something different with that of [56], where Nu 
is very sensitive to Pr (0.2 < Pr < 4). The fact that the power law of Nu is little sensitive to Pr 
as Pr > 4 in Fig. 2-8 & is consistent with the results of [55]. Here a simple explanation is presented for 
these results. 

Considering the steady state solution of Eq.(2-1) as Ra oo , we obtain the following equation 
by taking ^' = ^'Vk [48, 52]. 

P,-.J(V^*,*, = (_ + _) + R,_ p.2b) 

When Pr ^ 1 , then the left convection term Pr~^ J(V^^, in Eq.(2-26) can be ignored, and the 
buoyancy forcing term balances the viscous term V^\I' . In this case, the governing equations 

is independent of Prandtl number Pr , which is due to [52] . This is the reason why Nu and ^max are 
little sensitive to Pr for Pr > 10 in the numerical simulations by [55]. When Pr is order of 1 or even 
less than 1, the convection term can not be ignored and Pr plays a role in this case, as obtained by 
[56]. As Pr <C 1 , the buoyancy forcing should be balanced by strong convection at high Ra numbers, 
which can be seen from Eq.(2-16). This is also can explain the numerical results by [51], where the 
flow becomes unsteady and has strong eddy diffusion at Pr = 0.1 . 



§ 2.2 Partial-Penetrating Flow 

§ 2.2.1 Introduction 

Horizontal convection, in which the water is unevenly heated at the horizontal surface, was taken 
as a model of abyssal ocean circulation. The circulation, driven by density gradient and referred as 
thermohaline circulation (THC) at North Atlantic Ocean, is thought of an important energy con- 
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veyor belt and has great impact to climate change [42, 43, 44]. The horizontal convection become an 
important model to discuss the ocean energy balance [45, 46, 47]. Unlike the Rayleigh-Benard con- 
vection, the horizontal convection can be set to motion by any small temperature gradient. Moreover, 
the horizontal convection yields 1/5-power laws of Ra , comparing with the 1/4-power laws in the 
Rayleigh-Benard convection. 

The 1/5-power laws of Ra for flow strength (streamfunction maximum ^I'max ) and the heat flux 
(Nusselt number Nu), first found by [52], were later approved by both experiments [e.g. 52, 53, 54] 
and numerical simulations [e.g. 48, 55, 56, 66]. According to the scaling analysis [48, 52], there is 
a boundary layer near the surface, which is inverse proportion to 1/5 power of Ra . Both the flow 
strength and the heat flux are dominated by the scale of boundary-layer. Although this 1/5-power 
law of Ra is obtained for steady flow, it is still valid even for unsteady flow [e.g. 53] . 

The unsteady flow in horizontal convection was first found by numerically [51], then was observed 
in the experiment [53]. This unsteady flow is proved to be non-turbulent for that the energy dissipation 
turns to zero under the conditions of the viscosity v and thermal diffusivity k. being lowered to zero 
with Prandtl number Pr = zz/k being fixed [51]. Motivated by above investigations, [56] tried to find 
the bounded of Nu as Ra oo . However, their conclusion of Nu < cRa^^^ (for some constant c ) 
seems not tight enough as all the numerical simulations yield 1/5-power law. 

In a recent experiment, a new flow configuration referred as "parti-penetrating flow" was reported 
[54]. According to the measurement, the circulation cell is shallow and and no longer occupies the 
whole length of the tank. Though lots of numerical simulations on the horizontal convection, none 
of them have obtained such partial-penetrating flow. The main reason is that all of the former 
simulations used free slip condition on the walls [48, 55, 56, 62], the energetic circulation turns to be 
full-penetrating. While the laboratory experiments always require no slip on the walls [54], the viscous 
drag slows down the vigorous circulation. Noting that the power law fitted from the experimental data 
is somehow coarse, this can be improved by numerical simulations. Moreover, the partial-penetrating 
flow and the onset of such flow needs to reveal too. 

The main purpose of this paper is to investigate the partial-penetrating flow and to find a more 
accurate power law by numerical simulation, thus resulting in a more comprehensive view on this issue. 
The model and the scheme are sketched in § § 2.2.2, in which a benchmark solution is included. The 
establishment of the circulation and the partial-penetrating cell of horizontal convection are depicted 
in § § 2.2.3, with the onset of the partial-penetrating flow being discussed. Finally § § 2.2.4 concludes. 



§2.2.2 Model and Scheme 

We consider the the horizontal convection flows within the two-dimensional domain, and the 
Boussinesq approximation is assumed to these flows. The horizontal (y) and vertical (z) regimes are 
< y < i and < z < -D , respectively. Similar to [52] , the depth L is taken as reference length 
scale and A = D/L denotes the aspect ratio. We use A = 1 m present work as [52] did, while 
A = 0.675 is used in the experiments by [54]. Taking account of nondivergence of velocity field 
in Boussinesq approximation, the lagrangian streamfunction ^ and the corresponding vorticity to 
are introduced. The velocity 1? = (f , w) , where horizontal velocity v = ^ and vertical velocity 
u) = — ^ , respectively. The governing equations [48, 49, 50, 51, 56] in vorticity-streamfunction 
formulation are 
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(2-3a) 



Qy2 Q^2 

BT 

Vi{\/'^uj + Ra— ) 



(2-3b) 



(2-3c) 



where J {^,4)) = ~ denotes the nonhnear advection term. There are two important di- 

mensionless parameter in Eq.(2-3), i.e. Rayleigh number Ra = UT^TgL^ / {kv) and Prandtl number 
Pr = i^/k , where g , ax , AT , L , k and u are gravity acceleration, thermal expansion coefficient, 
surface temperature difference, length of horizontal domain, thermal diffusivity and kinematic vis- 
cosity, respectively. The boundary condition is the same with the experiment: the surface buoyancy 
forcing is T = sin(|y) , and no slip boundary condition is applied to walls except for surface. 

There are two important quantity describing the circulation, i.e. the non-dimensional stream- 
function maximum and the non-dimensional heat flux. The non-dimensional streamfunction maxi- 
mum ^'max = ^max/^ i where ^'max the maximum of the dimensional streamfunction. The non- 
dimensional heat flux is defined as Jt = dT/dz on the heated surface. Nusselt number Nu , which 
is defined here as the maximum of dT/dz on the top surface. This definition of Nu is something 
different from the others [e.g. 53, 56]. 

The above Eq.(2-3) is solved with finite different method in non-uniform grids. Crank-Nicholson 
scheme and Arakawa scheme [e.g. 60, 61] are applied to discretize the linear and nonlinear terms, 
respectively. Comparing to the other schemes, Arakawa scheme is more accuracy but more expensive, 
and it has also been applied to horizontal convection flows at high Rayleigh number [62]. 

To validate the scheme, we calculate the nature convection problem with the resolution of 80 x 80 
meshes. Table 2-1 compares the present solutions with bench mark solutions from [63], who solved 
the problems using a pseudo-spectral Chebyshev algorithm and [64], who used fourth-order compact 
finite difference scheme. It is obvious that the our scheme is very accuracy for such fiows. 

§2.2.3 Results 

1. spatial resolution 

First, we investigate grid dependency of the solutions to ensure that the numerical simulations 
are valid. The boundary condition is the same with the experiment: the surface buoyancy forcing is 
T = sin(^y) , and no slip boundary condition is applied to walls except for surface. To this purpose, 
a case of Ra = 2 x 10^ and Pr = 1 is calculated with grids of three different resolution, i.e. the 
horizontal number of meshes = 40 , A = 64 and A^ = 80 . We find that the resolution of grids 
must be fine enough, otherwise some unphysical time-depend solutions would be obtained. 

Fig. 2-9 a depicts the time evolution of the maximum ^'max • The solutions tend to be steady as 
time t > 1 for A = 64 and A = 80 . While it becomes time-dependent for A = 40 . It implies 
that some unphysical time-dependent solutions might be obtained if the spatial resolution is not fine 
enough. To exclude the unphysical time-dependent solutions, the numerical simulations must be 
obtained with sufficient spatial resolution which depends on the Rayleigh number Ra . As Fig. 2-96 
shows, the minimal number of horizontal meshes A^ is to obtain correct results directly proportion to 
Ra^''^ . Taking account of Ra oc , this means N (x L: the longer L is, the larger N is. To obtain 
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the physical solutions, must be within the stable regime in Fig. 2-9 6. According to our calculations, 
the flow is still steady and stable for Ra < 10^^ . 

It is from Fig.2-96 that Ay = L/N = CRiKi^y/^ /{axATgy/^ , where Ay and Cr = 10 are the 
mesh size in y direction and the coefficient, respectively. The smaller n and v are, the smaller the 
mesh should be. For the molecular kinematic viscosity v = 1.5 x 10~^ cm^/s and thermal diffusivity 
K = 1.3 X 10~^ cm^/s in the case of "run 16" by [54], the mesh Ay should be 2.1mm, which is 
smaller than Kolmogorov scale r/ = (z/'^/e)^/^ = 5.8 mm , where e = 2 x 10~^cm^/s is dissipation rate 
in the field [54]. So this implies that the mesh should be fine enough to resolute Kolmogorov scale 
eddies. 

The resolution requirement of Ay implies that the bound of L in numerical simulations is about 
laboratory scale if molecular viscosity and diffusivity are used. And eddy viscosity and diffusivity are 
required, when the length of ocean scale is considered. 



2. the establishment of circulation 

Then the circulations are obtained numerically for Ra > 10'' . Similar to [51], the dimensionless 
time t = t* /t , where t* and r = D"^ / k, are the dimensional and the unit scaling times, respectively. 
A typical value of i = 1 is about 80 hours in the dimensional time, given D = 20 cm and n = 
1.4 X 10^'^ cm^/s , which are approximate to the values used in the experiments by [53, 54]. In the 
following simulations, we use the total kinetic energy and ^max as indexes to oversee the evolution 
of the circulations, where Ej^ = §{y'^ + w'^)/2dydz. It is found that the flow is established very fast, 
but it takes a very long time for the flow to approach to equilibrium state. 

To illuminate this, the circulation of Pr = 8 and Ra = 5 x 10^ is taken as an example. Fig.2- 
7 displays the time evolutions of total kinetic energy E^ and ^'max with time t . The circulation 
is established very fast, and both Ef^ and ^'max reach their 95% of the equilibrium values within 
t = 0.05 (4 hours). In addition. Fig. 2-76 agrees well with the experiments (e.g. Fig. 2 in [54]), and 
this fast establishment of the circulation was also noted by [53]. However the total time to get the 
quasi-equilibrium state is relatively longer, it takes at least te = 2.5 (200 hours) for the flow being 
steady in present simulation. During this process, the flow near the surface is fully established and 

(a) (b) 




I N 

Figure 2-9: (a) The maximum of streamfunction ^max vs time t for Ra = 2 x 10^ . The solid, dashed 
and dash-doted curves are solutions with = 40 , = 64 and = 80 , respectively. (6) The stable 
and unstable regime on the plot of Rayleigh number Ra vs A^ . 



34 



§2.2 Partial-Penetrating Flow 



(a) 



0.8 - 



0.6 




0.4 



0.2 - 



IT - - -0.005 _ 

0- 



-0- 



-0.0002 - 

_1j I I 



0.2 0.4 0.6 0.8 1 

3' 




0.6 



0.4 



0.2 



(b) 



-0.375' 



-0.01- 





0.2 0.4 0.6 0.8 1 

y 





Figure 2-10: The flow fields (streamfunction at four different time steps: t = 0.5 (a), t = 1.0 (6), 
t = 1.5 (c) and t = 2.5 {d) at Pr=8 and Ra = 5 x 10^ , solid curves for clockwise flow ( ^' > ) and 
dashed curves for anticlockwise flow < 0), respectively. The partial-penetrating cells > 20) 
are shadowed and the counter intervals are 10 for ^' > in each figures. 
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approaches to equihbrium state within t = 1 , as Fig.2-10a and Fig. 2-106 show. Meanwhile, there are 
several very weakly clockwise and anticlockwise circulation cells below the primary circulation cell, 
which are secondary flows due to the drag by upper primary circulation. These secondary cells were 
also observed in the experiments by [54]. As the time goes on, the primary circulation cell becomes 
stronger and larger, and the secondary circulations are weaker and smaller. So that the counter line 
of ^' = 10 becomes deeper and deeper as shown from Fig.2-10a to Fig.2-10d Finally, the primary 
cell fulfills the whole tank at t = 2.5 (Fig.2-10(i). 

Then, the flows within 10^ < Ra < lO^'' are calculated. Fig.2-8a shows the 1/5-power law of 
Ra for ^'max at different Prandtl numbers (e.g. Pr = 1 , Pr = 4 , Pr = 8 and Pr = 10). Noting 
that Rup = RaPr^^'^ instead of Ra are used in Fig.2-8a, it implies that the flow is dominated by 
thermal diffusivity n . The larger Prandtl number is, the stronger the flow is. All the direct numeric 
simulation (DNS) data at Pr = 1 , Pr = 4 , Pr = 8 and Pr = 10 yield, 

^'max = 0.75i?ay^ = 0.75i?ai/^Pr^/^° (2-4) 

On the other hand, [54] fitted their experiments data into two different 1/5-power laws: bigger one 
for the reduced acceleration due to gravity g' > 2.5 cm/ s'^ , the smaller one for g' < 2.b cm/ . 
Equation (2-4) is similar to but more accurate than the bigger one obtained by [54]. It is notable 
that all the experimental data lie around this line (right triangles in Fig.2-8a), as the Prandtl number 
is considered. Taking account that the aspect ratio in DNS is not exact the same as that in their 
experiments, these data are consistent well with each others. Fig. 2-86 shows the 1/5-power law for 
Nu . This is consistent with the scaling analysis, i.e. the thermal boundary-layer is inverse proportion 
to Ra^/^ . It is obvious that Nu is less sensitive to Pr than that of ^'max i so that Nu seldom 
changes for Pr > 4 . 



3. the partial-penetrating cell 

It is notable that the main circulation near the surface is seldom changed during the establishment 
process, especially the close circulation cell shadowed in Fig. 2-10. As this cell is shallow and the 
circulation is only near the surface, it is referred as the "partial-penetrating cell" after [54]. For 
example, the shadowed cell height is about 1/6 of the total depth in Fig. 2-10, so that the whole cell 
is within the boundary-layer near the surface. This kind of flow also exist in other Rayleigh numbers, 
where the shallow and close cells like that in Fig. 2-10. 

An objective definition of "partial-penetrating cell" is more convenient for further discussion. In 
present investigation, the penetrating depth of the cell Dc is defined as the depth of a close circulation 
cell, which contains 60% of the total amount. Here the concept of partial-penetrating cell has such 
meanings: (1) above all, the close circulation cell is very shallow comparing to its width, so that Dc 
is within the boundary-layer near the surface, (2) consequently, the close cell is seldom affected by the 
bottom boundary, (3) moreover, the flow in close circulation cell is dominant of the main circulation, 
e.g., the flow rate in the shallow cell is about 60% of the total amount (e.g. = 20 in Fig. 2-10). 
On the contrary, "the full-penetrating circulation" is referred to the circulation takes whole depth of 
the tank hereinafter. 

Using the above definition, the shadowed cell in Fig. 2-10 is the partial-penetrating cell. First, 
it is obvious that the shadowed cell satisfies conditions (1) and (3). Second, the close cell is seldom 
affected by the bottom boundary. To illuminate this, we descend the depth of the water tank (or 
equally descending aspect radio ^4). Fig. 2-11 clear depicts that the partial-penetrating part is little 
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Figure 2-11: The flow fields (streamfunction ^ ) near the forcing surface of three respective aspect 
ratios: A = 1 (a), A = 0.5 (6) and A = 0.25 (c) at Pr=8 and Ra = 5x 10* . The partial-penetrating 
cehs are shadowed and the counter intervals are 10 in each figures. 

sensitive to bottom, if the depth of the tank is large enough {Dc < D < 2Dc in Fig.2-llc). It may 
also be noted that the flow in Fig.2-llc is very similar to the flow field near the surface ( 0.75 < z < 1) 
in Fig.2-10a. For that at that time in Fig.2-10a, the boundary of main circulation 4' = is about 
z = 0.75 , hence the flow near the top is approximation of the case in Fig.2-llc. This also implies that 
the shadowed cell in Fig. 2-10 is partial-penetrating. 

As the partial-penetrating cell is near the forcing surface, it is quite independent of the flows 
in the middle and the bottom. First, the partial-penetrating cell changes little even when the full- 
penetrating circulation emerges in Fig.2-10c, d. Second, In addition, Fig.2-11 clear depicts that 
the partial-penetrating part is little sensitive to bottom, if the depth is large enough. Thus the 
partial-penetrating cell ( ^' > ^'c ) is unchanged given the tank is deep enough. On the contrary, the 
full-penetrating part is very sensitive to the bottom, as the flow this part of flow approaches to the 
bottom. 

4. the constrain for partial-penetrating flow 

However, not all the flows have partial-penetrating cells like that in Fig. 2-10, e.g. the flow at 
Pr=l and Ra = 5 x 10* in Fig. 2-12. It is obvious that the circulation of Pr = 1 , which is mainly 
full-penetrating, is more deep comparing to Pr = 8 . So only the close circulation cell for ^' > 35 is 
shallow enough, and it is about 1/6 of the total depth (Fig. 2-12). While this shallow cell is neither 
dominant in the whole circulation nor immune to the changes of bottom boundary, which is shown 
in Fig. 2-13. Even this small shallow cell is sensitive to the bottom boundary. Compare this with 
Fig.2-11, the close circulation cell for ^' > 35 is not like the partial-penetrating cell. 
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Figure 2-12: The flow fields (streamfunction ^' ) for full-penetrating flow at Pr=l and Ra = 5 x 10^ . 
The counter interval is 5, and ^ > 35 is shadowed. 



(a) 




Figure 2-13: The flow fields (streamfunction ^ ) near the forcing surface of three respective aspect 
ratios: A = 1 (a), A = 0.5 {b) and A = 0.25 (c) at Pr=l and Ra = 5 x 10*^ . The streamfunction 
^' > 35 are shadowed and the counter intervals are 5 in each figures. 
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Figure 2-14: Dc vs. Ra. The solid, dashed and dash dotted lines are power laws of Ra respectively 
for Dc at Pr = 6 , Pr = 4 and Pr = 2 . 

Then a new problem emerges. As there might be full-penetrating flows, what's the constrain for 
the partial-penetrating flow? To investigate this, Dc is employed as an index of partial-penetrating 
flow according to the above definition. If Dc satisfies —1/5 power law of Ra , then the flow is 
within the boundary-layer and the flow is partial-penetrating, for that the boundary-layer satisfies 
— 1/5 power law of Ra [48, 52, 53]. Else if Dc does not satisfy —1/5 power law of Ra , then the 
flow is full-penetrating. Thus, according to the power law for Dc , the flow can be referred as either 
partial-penetrating or full-penetrating. 

In fact, the critical parameter governing the exitance of partial-penetrating cell is Prandtl number 
Pr . The larger the Prandtl number is, the more obvious the partial-penetrating cell is. It is from 
Fig. 2-14 that the smaller Pr is, the bigger or deeper Dc is. Only for the flows at Pr > 6 , does Dc 
satisfy the —1/5 power law of Ra . Dc seldom changes with Ra even for the flows at Pr = 1 . 

It is notable that partial-penetrating flows exist for Pr > 6 , this can be understood from scale 
analysis [48, 52]. When Pr » 1 , then the left convection term Pr^^ J(V^^', in Eq.(2-36) can 
be ignored, and the buoyancy forcing term Ra^ balances the viscous term . In this case, the 

governing equations is independent of Pr and the boundary layer is inverse proportion to 

When Pr is order of 1 or even less than 1, the convection term can not be ignored. So that the strong 

convection is full-penetrating as depicted in Fig. 2-12. 

In words, it is Pr that governs the flow pattern. When Pr > 6 , the flow is partial-penetrating. 
While Pr < 4 , the flow is full-penetrating. 

§ 2.2.4 Conclusion 

The experiments of horizontal convection in a rectangle cavity are simulated numerically. The 
simulations agree well with the experimental data, and a more extensive 1/5-power law of Ra is 
obtained by fitting the DNS data. The partial-penetrating flows are revisited by numerical simulations. 
It is Pr that governs the the existent of the partial-penetrating flow. When Pr > 6 , the flow is partial- 
penetrating. While Pr < 4 , the flow is full-penetrating. 
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§2.3 Power Laws 

§ 2.3.1 Introduction 

The power laws in horizontal convection is one of the fundamental for understanding the flows. 
It is reported from [52] and [54] that the flow strength (the maximum of stream function ^max ) is 
always increasing with Ra}/^ . However, the flux of the flow can't be measured very accuracy, and 
the parameters are bounded within Ra < 10^'^ due to experimental inconvenientness. Moreover, the 
numerical simulations of [51] ( 0.1 < Pr < 10 ), [56] ( 0.2 < Pr < 4 ) can obtain both the flow strength 
and the heat Hux. But the Rayleigh number Ra used in numerical simulation are lower than the 
maximum value used in the experiments. It is found that ^max ~ Ra^^^ to ^max ~ Ra}^"^ , and 
Nu oc Ra}^^ [56], which is consist with the results by [62] at Pr = 10. The flow strength depends 
continuously to Ra from ^max oc Ra to ^max oc Ra^/^ . 

On the other hand, Paparella and Young (2002) [51] have claimed that the Hux trends to zero 
as thermal diffusion trends to zero and Pr is fixed. Besides, SIggers et. al (2004) obtained that the 
maximum of heat flux Nu is proper to Ra^/^ , which is not supported by their numerical simulations. 
Then two possibilities are suggested. First, their theory is not accuracy. Second, the numerical 
simulation is not accuracy. 

Motivated by the above problem, the power laws of horizontal convection is investigated. Similar 
to the previous investigations, we consider the horizontal convection flows within the two-dimensional 
domain with aspect ratio A = 0.1 . In the following parts, the flows at very high Rayleigh numbers 
are simulated, and the power laws are obtained by fitting the individual result. 

§ 2.3.2 Results 

When the Rayleigh number is relatively lower, the flow is weak and the heat transportation is 
dominated by heat conduction [48, 56]. As the Rayleigh number increases, the thermal convection 
enhances to balance the surface heat flux. Fig. 2- 15 shows the stream function and temperature field 
of Ra = 10^ . There is a weak upwelling How in the center and two strong downwelling flow in the 
two sidewalls, which are driven by central heating and sidewall cooling. Accordingly, the temperature 
field deforms it horizontal isothermal lines as the convection. There is a warm tongue (shadowed in 




Figure 2-15: The flow stream function (a) and temperature field (b) of Ra = lO'^ . It is steady and 
stable and symmetric with middle plume forcing, solid and dashed curves for positive and negative 
values, respectively. 
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Fig.2-15b), which is due to the vigorous central flow near the top surface and the central flow near the 
bottom bringing the cold water to the middle. For example, the isothermal line of T = 0.3 extends 
to cold water near top but but to warm water near bottom. 

When the Rayleigh number is higher, the flow field is more complex. Fig. 2-16 depicts the flow 
filed and temperature field of Ra = 10^*^ , which are something different from the flow of Ra = 10^ . 
First, the downwelling flow departs into two separate parts: one descends the middle, the other can 
descend to the bottom. Thus there are two vortex centers in the main circulation as in Fig.2-16a. 
Besides, there is a reverse secondary circulation along the bottom near the sidewalls (1/10 departs 
from sidewall), which have no mass interchange with the main circulations (dead water). Consequently, 
they have no contribution to the main thermal convection. Similar to Fig.2-15a, the isothermal lines 
are concentrated and horizontal in the center. Moreover, the isothermal lines in Fig.2-16b shows 
there are a temperature boundary layer and a thermocline near the surface. In total, the flow field 
are dominated by cold water, where T < 0.2 except for the boundary layer. This contracts to the 
Fig.2-15b. 

Compare the temperature fields in Fig.2-15b and in Fig.2-16b, the temperature boundary emerges 
near the top surface as the increase of Ra. The temperature is homogenous in the middle and 
bottom, except for near the sidewalls. The stable stratification also emerges during this process. 
It is horizontal temperature gradient that driven the surface circulation, which can be seen from the 
flow flled. However, the bottom circulations are dominated by kinetic adjustment, which has nothing 
to do with temperature distribution. 



§2.3.3 Power laws 

The above investigation have obtained the flow pattern of two different Rayleigh numbers. The 
general properties of the horizontal flows are addressed here, especially for the sensitivity of Ra. To 
this purpose, the power laws for flow strength and heat flux are calculated. 

First, the flow strength vs. Ra is investigated, which is the main parameter of intensity in the 
thermohaline circulation. It is from Fig.2-17a that the flow strength ^max is linear proper to Ra as 
Ra < 10^ , and that ^max is 1/5-power proper to Ra as Ra > 10^ , between is the transition regime. 




Figure 2-16: The flow stream function (a) and temperature field (b) of Ra = 10^^ . It is steady and 
stable and symmetric with middle plume forcing, solid and dashed curves for positive and negative 
values, respectively. 
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The power law of Ra at high Rayleigh number is, 

^max = 0.Q2Ra^/\ (2-5) 

Similar to the numerical simulation by Siggers et al (2004), "^max ~ Ra^^^ to ^max ~ Ra^/^ as 
10^ <Ra< 10^ . But there is an obvious 1/5-power law as Rayleigh number is higher enough, which 
implies that this 1/5-power law is valid only at relatively higher Ra. The power law obtained by 
Siggers et al (2004) is only valid for transition regime. 

Second, the power law for heat flux is also studied. As in Fig.2-17b, the heat flux is seldom 
changed as the increase of Ra. In fact, the heat flux is dominated by heat conduction in this regime 
[48, 56]. However, this changes to convection dominance as Ra larger than a critical value Ra ~ 10^ , 
and Nu increases quickly as Ra . Then, a 1/5-power law emerges as Ra > 10^ , 

Nu = 0.75i?a^/^ (2-6) 

Similar to the numerical simulation by Siggers et al (2004), there is no 1/3 -power law regime 
in the simulations, though Ra is much larger here. The simulations supports their first hypothesis 
that the theory is overestimation to the heat flux. According to the analysis by Rossby [52], and 
Quon [48], the thickness of temperature boundary layer is proper to 1/5-power of Ra, which implies 
Nu ~ Ra^ /5 . Compare the power laws for ^max and Nu , there are always three regimes: linear 
regime, transition regime and 1/5-power regime. However, the transition from one regime to an other 
is not synchronous for ^max and Nu , where the transition for heat flux is leading than for the flow 
streamfunction. 



§ 2.3.4 Summary 

In summary, the horizontal convection at high Rayleigh number in a rectangle cavity with aspect 
ratio of 1 : 10 is numerically simulated. According to the results within the regime of 10^ < Ra < 
10^^ , three continues regimes are obtained: linear regime (10^ < Ra < 10^), transition regime 
( 10^ < Ra < 10^ ) and 1/5-power law regime ( 10^ < Ra < 10^^ ). For the flow strength, a 1/3-power 




Figure 2-17: The flow stream function ^max (a) and heat flux (b) vs. Ra. 



42 



S 2.4 Onset of Instability 



law of Ra is fitted when Ra is not high enough ( 10'' < Ra < 10^ ). However, a 1/5-power law is 
obtained as Ra is high enough ( 10^ < Ra < 10^^ ). The 1/5-power law confirms Rossby's analysis 
and implies that 1/3-power law of Ra for Nusselt number by Siggers et al. is over estimation. 

§ 2.4 Onset of Instability 

§ 2.4.1 Introduction 

Horizontal convection, in which the water is unevenly heated at the horizontal surface, was taken 
as a model of abyssal ocean circulation. As the abyssal ocean circulation plays an important role in 
climate change, the horizontal convection has intensively been explored in recent years [51, 53, 54]. It 
can be set to motion by any small temperature gradient, unlike the Rayleigh-Benard convection. But 
similar to Rayleigh-Benard convection, the horizontal convection may be unsteady at high Rayleigh 
numbers Ra . There is a critical Rayleigh number Rac , and the steady flow is unstable and becomes 
unsteady when Ra > Rac ■ The unsteady flow in horizontal convection was first found by numerical 
simulation [51], then was observed in the experiment at Ra > 10^^ [53]. This unsteady flow is proved 
to be non-turbulent even as Ra — > oo , though the flow field seems to be chaotic [51]. The investigation 
on the unsteady horizontal convection flow is relatively less, except for [51, 53, 67]. However, they 
have mainly focused on how the turbulent plume maintains a stable stratified circulation. Yet how 
the horizontal convection turned to be unsteady remains an elusive problem. 

To understand this problem, both Rac for the onset of unsteady flow and instability mechanism 
are of vital. Paparella and Young [51] found Rac ~ 2 x 10® at Pr = 1 in their simulations, which is 
significantly smaller than others' results. For example, Rossby (1965), Wang and Huang (2005) found 
the flow is steady and stable for Ra < 5 x 10® in their experiments [52, 54]. Yet some other numerical 
simulations [55, 56, 62] have not found unsteady flows for Ra < 10^ . Paparella and Young [51] 
explained this difference as: (i) lower aspect ratio ( H/L = 1/4 ) than the experiments and (ii) middle 
plume forcing instead of sidewall plume forcing in the experiments. Both may lead to destabilization of 
the flow at lower Rayleigh numbers. However, their hypotheses have not been intensely investigated. 
According to a recent investigation, the flow is still stable for Ra < 10^^ even at a much lower aspect 
ratio {H/L = 1/10) [68]. Thus, it maybe the middle plume forcing that leads to destabilization at 
lower Rayleigh numbers. 

Our interest here is to verify their second hypotheses. Is the flow with middle plume forcing less 
stable than the sidewall plume forcing? How the instability occurs? To investigate these problems, 
more accurate numerical prediction of Rac is need for both forcing cases, for the spatial resolution of 
simulation is very coarse used (e.g. 128 x 32 meshes are used in [51]). Then the flow field under both 
middle and sidewall plume forcings are compared, which leads to an affirmative answer of the above 
problem. 

Similar to the previous investigations, we consider the horizontal convection flows within the 
two-dimensional domain, and the Boussinesq approximation is assumed to be valid for these flows. 
As shown in Fig. 2-18, the horizontal (y) and vertical (z) regimes are < y < L and < z < H , 
respectively. Similar to [52], the depth L is taken as reference length scale and A = H/L = 1/4 de- 
notes the aspect ratio. Taking account of nondivergence of velocity field in Boussinesq approximation, 
the Lagrangian streamfunction ^ and the corresponding vorticity uj are introduced. The velocity 
"u = {v,w) , where horizontal velocity v = ^ and vertical velocity w = , respectively. The 
governing equations in vorticity-streamfunction formulation are [48, 51, 56]: 
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dT , , ,d'^T d'^T^ 

^ + J{^, io)=- FriV^LO + Ra— ) (2-7b) 
V^^r = -co (2-7c) 

where J{^,(f)) = — denotes the nonhnear advection term. There are two important di- 

mensionless parameter in Eq.(2-7), i.e. Rayleigh number Ra = ax^TgL^ /{nv) and Prandtl number 
Pr = , where g, ot , AT, L, k and v are gravity acceleration, thermal expansion coefficient, 
surface temperature difference, length of horizontal domain, thermal diffusivity and kinematic vis- 
cosity, respectively. Alternatively, Paparella and Youngs used vertical length H as length scale, so 
Ra = 64:RaH , where Ran is the vertical Rayleigh number by using vertical length H as unit [51]. 

More specifically, we consider the horizontal convection in a rectangle tank at Pr = 1 . The tank 
has same velocity boundary condition as that in [51], i.e. free slip and no stress at the walls. In 
addition, two different surface forcings are used, which are central symmetric. One is middle plume 
forcing as T = [1 + cos(27ry)]/2 [51], the other is sidewall plume forcing as T = [1 — cos(27ry)]/2 
[48, 68]. Comparing these with one cell forcing T = cos(7ry/2) [55], there are two symmetric cells in 
the flow field under such forcings (e.g. Fig. 2-18 and Fig. 2-24), when the flow is symmetrically steady 
and stable. In addition, the middle plume forcing in the left cell is the same with the sidewall plume 
forcing in the right cell (see Fig. 2-18 behind). Thus in the steady flows, both forcings will lead to the 
same flow patterns except for a position shift, which is proved by the following investigation. 

There are two important quantity describing the circulation, i.e. the non-dimensional streamfunc- 
tion maximum and the non-dimensional heat flux. The non-dimensional streamfunction maximum 
^max = ^ma,x/^ ) where ^'max is the maximum of the dimensional streamfunction. 

The above Eq.(2-7) is solved with finite different method in non- uniform grids. Crank-Nicholson 
scheme and Arakawa scheme [e.g. 60, 61] are applied to discretize the linear and nonlinear terms, 
respectively. Comparing to the other schemes, Arakawa scheme is more accuract but more expensive, 
and it has also been applied to horizontal convection flows at high Rayleigh number [62, 68]. Table 
2-1 shows the validation of the scheme with nature convection problem. A fine spatial resolution mesh 
of 512 X 128 is used to eliminate numerical instability. 



§ 2.4.2 Results 

First, the middle plume forcing is considered, which is steady and stable for Ra < 5.5 x 10^ . 
Fig. 2-18 shows the flow field (a) and temperature field (b) of Ra = 5 x 10^ with ^max = 59.83 , in 
which the flow is symmetric, steady and stable. In this case, the center line symmetrically separates 
the flow field into two parts, like a free slip wall. There is a vigorous downward jet in the center of 
tank corresponding to the middle plume forcing (Fig.2-18a), where the vertical velocity field has a 
minimum of w = —2513 (Fig.2-19b). The center jet leads to the clockwise and anticlockwise plume 
cells in the left and right part of tank, respectively. In the left circulation cell, the flow sinks quickly 
along the center line and upwells clockwise along the left side wall with relatively slower speed, which 
can be also seen from the vertical velocity of the flow (Fig. 2-19). Besides, there are two horizontal 
jets respectively near top and bottom walls in the left circulation cell (Fig.2-19a). Totally, there are 
2 horizontal jets near wall and a vertical jet at the center in each cell. Contract to the flow field, the 
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Figure 2-18: The flow stream function (a) and temperature field (b) of Ka = 5 x 10^ . It is steady and 
stable and symmetric with middle plume forcing, solid and dashed curves for positive and negative 
values, respectively. 

temperature field is very simple. An obvious boundary layer exists near the surface in temperature 
field, which leads to a 1/5-power law of Ra for heat flux [e.g. 48, 52, 56]. And below the temperature 
boundary layer, the temperature is almost homogeneous due to the convection. Thus there is a very 
strong stratiflcation near the surface ( dTjdz ~ Ra^l^ ) but a very weak stratiflcation in other region 
{dT/dz ~ 0). As the above case is stable, so the critical Rayleigh number must be larger than 
5 X 10^ , which is signiflcantly larger than the value obtained before [51]. 

To flnd the critical Rayleigh number Rac , the growth rate of perturbation 4){t) is calculated 
numerically. And (/>(i) is assumed to satisfy 4){t) = 6^^*0(0) , where a = ar+itJi is the complex growth 
rate of disturbance. It is found that the onset of unsteady flow is at Rqc = 5.5377 x 10^ , as shown in 
Fig. 2-20. For Ra = 5.53 x 10^ , the flow is stable and the growth rate is approximately ar = —0.12 . 
But the flow is unstable and the growth rate is approximately ar = 0.03 for Ra = 5.54 x 10* . Thus 
the critical Rayleigh number Rac is obtained 5.53 x 10* < Rac < 5.54 x 10* . The accurate value 
of Rac = 5.377 x 10* is obtained by interpolating from the above result. Moreover, the onset of 
unsteady flow is found to occur via Hopf bifurcation. As Fig. 2-20 shows, the image part of growth 



(a) 









■"::v:^(®CE:r. --' 
















"—"'■'''^-506 






-1000 
















Figure 2-19: The horizontal (a) and vertical (b) velocity flelds of Ra = 5 x 10* . It is steady and 
stable and symmetric with middle plume forcing, solid and dashed curves for positive and negative 
values, respectively. 
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rate is nonzero and the eigenmode of perturbation is periodic. This Hopf bifurcation of the horizontal 
convection has not seen reported yet, and privous investigations dealt only with chaotic flows. 

Meanwhile, the evolution of the perturbational vorticity fields during the fist half period at t = 
(a), t = T/8 (b), t = r/4 (c) and t = 3T/8 of Ra = 5.54 x 10^ are depicted in Fig.2-21, respectively 
The perturbational vorticity fields are symmetric about centerline, which implies that the horizontal 
velocity is nonzero at centerline. It can be seen that the perturbation tripole A (the shadowed ellipse 
in Fig.2-21a) is generated from central downward jet, then propagates and amplifies along the central 
jet downward to the bottom wall (Fig.2-21b,c,d). When tripole A approaching to the bottom, it 
becomes weaker and weaker and breaks into two parts: the left and the right near the bottom, which 
can be seen from the evolution of tripole B (the shadowed rectangle in Fig.2-21a). And the mean flow 
advects the broken vortexs horizontally along the bottom wall (Fig.2-21b,c,d). Then in the second 
half period, a reverse tripole will generate right the same place of vortexes A at t = T/2 , and the 
same story repeats for it, which is shown in Fig. 2-22. In short, the perturbations generate and amplify 
in the central vertical jet, but are propagated and weakened along the horizontal wall. 

Further investigation shows that the instability of flow occurs due to shear. First, as we noted, the 
instabilities always occur in the center and propagate along the mean flow. Second, this trigger place 
locates in the area where there is a vigorous jet with strong shear (see e.g. Fig.2-19b and Fig. 2-23) 
(Fig. 2-23). As the stratification is very weak here (see e.g. Fig.2-18b), so that the flow in this region 
is dominated by momentum dynamics other than thermal dynamics. All these imply that the onset 
of instability leading to unsteady flow is due to shear instability at larger Raleigh numbers, which is 
much different from Rayleigh-Benard instability. However, shear is not the only sufficient condition for 
instability. For example, the instability near the top surface is suppressed due to strong stratification 
(Fig.2-18b), though both the velocity (Fig.2-19a) and the shear (Fig. 2-23) near top surface are still 
very large. In words, the onset of instability is due to velocity shear (shear instability) other than 
thermally dynamics (thermal instability). 

Second, the sidewall plume forcing is considered. Fig. 2-24 shows the flow field and temperature 
field of Ra = 5 x 10^ , in which the flow is symmetric, steady and stable like that under middle plume 
forcing. There are two strong downward jets near the walls corresponding to the sidewall plume 
forcing (Fig.2-24a). As mentioned above, the sidewall plume forcing will lead to exactly the same flow 
pattern as the middle plume forcing does except for a position shift, which can be seen from Fig. 2-18 
and Fig. 2-24. As the flow is stable, the center line like a free slip wall symmetrically separates the 
two cells. The left cell in Fig. 2-18 is exactly the same with the right cell in Fig. 2-24. However, the 




Figure 2-20: Growth rate 0"^ (solid) and cjj (dashed) vs. Ra , respectively. 
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Figure 2-21: The perturbational vorticity fields at t = (a), t = T/8 (b), t = T/A (c) and t = 3T/8 
of Ra = 5.54 x 10^ , solid and dashed curves for positive and negative values, respectively. 
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Figure 2-22: The perturbational vorticity fields at t = T/2 (a), t = 5T/8 (b), t = 6T/8 (c) and 
t = 7T/8 of Ra = 5.54 x 10^ , solid and dashed curves for positive and negative values, respectively. 
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Figure 2-23: The vorticity of Ra = 5.54 x 10^ with vertical velocity w (shadowed as w > 1800), 
solid and dashed curves for positive and negative values, respectively. 

flow is much more stable with the sidewall plume forcing. And the critical Rayleigh number is found 
Rac ~ 1.85 X 10^° (with 768* 192 meshes) in this case. As noted above, the flow is much more stable 
with the sidewall plume forcing than that with the middle plume forcing, though both forcings lead 
to the same flow patterns. This is very interesting, and can be understood from the mechanism of 
instability. 

It's found that the rigid wall suppresses the perturbation, which leads a more stable flow with the 
sidewall plume forcing than that with the middle plume forcing. As the flow loss stability is due to 
strong velocity shear in the center in horizontal convection, the smaller of the shear the more stable of 
the flow. In the case of middle plume forcing, the perturbation with nonzero horizontal velocity occurs 
at strongest downward jet. And the perturbed flows cross the center line and propagate downstream. 
However, in the case of sidewall plume forcing, these crossing flows are suppressed by rigid walls. So 
that the critical Rayleigh number is much larger in this case. Paparella and Young (2002) hypothesized 
that middle plume forcing may lead to a destabilization of the flow. Here this hypotheses is proved 
both physically and numerically. 



§ 2.4.3 Conclusion 

In conclusion, the onset of unsteady flow is found to occur via a Hopf bifurcation in the regime 
of Ra > Rqc = 5.5377 x 10^ for the middle plume forcing at Pr = 1 , which is much larger than 
the previously obtained value. Besides, the onset of unsteady flow is due to shear instability of 
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Figure 2-24: The flow field (a) and temperature field (b) of Ra = 5 x 10® , which are steady and stable 
and symmetric with sidewall plume forcing. 
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central downward jet. Finally, the second hypotheses of Paparella and Young (2002) for instability 
is numerically approved, i.e. the middle plume forcing can lead to a destabilization of the flow at 
relatively lower Rayleigh numbers. 
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